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1.  INTRODUCTION 


During  the  past  five  years  the  microcomputer  has  developed  into 
a  powerful  computational  aid  to  the  engineer.  Currently,  Z80-based 
microcomputers  are  capable  of  performing  most  day-to-day  numerical- 
fluid-mechanics  tasks  such  as  data  preparation/reduction,  two-di¬ 
mensional  boundary- layer  computations  with  very  fine  finite-differ¬ 
ence  grids,  and  many  other  practical  computations  including,  on 
relatively  crude  finite-difference  grids,  time-averaged  Navier-Stokes 
solutions.  With  the  exception  of  the  latter  example,  these  compu¬ 
tations  usually  can  be  done  as  fast  on  a  microcomputer  as  running  on 
a  scientific  mainframe  computer,  if  you  include  the  time  required 
for  transmitting  the  input  data  and  final  printout  via  a  terminal 
Interface. 


This  is  true  because  of  the  overall  real  time  spent  using  a  main¬ 
frame  computer  with  the  microcomputer  as  a  terminal,  typically  only 
52  is  devoted  to  the  actual  execution  of  the  program.  The  other 
952  consists  of  signing  on,  transmitting  input  data,  and  eventually 
transmitting  the  final  printout.  Obviously  the  time  distribution 
varies  depending  on  the  overall  execution  time  and  the  amount  of 
print  taken.  However,  our  experience  over  a  two-year  period  has 
shown  a  52-952  split  to  be  quite  representative  for  numerical-fluid- 
mechanics  applications. 


On  the  basis  of  the  great  power  of  modern  microcomputers,  this  pro¬ 
ject  has  been  undertaken  with  the  overall  aim  of  developing  a  micro- 
computer-based  computational  procedure  for  predicting  viscous  flow 
about  a  ship  hull.  The  theoretical  foundation  of  the  procedure  is 
(a)  classical  thin-ship  theory^"  for  the  inviscid  velocity  distri¬ 
bution  on  the  ship  hull  and  (b)  a  three-dimensional  boundary-layer 

p 

program,  EDDY3,  which  incorporates  the  Wilcox-Rubesin  two-equation 
turbulence  model.  In  a  straightforward  manner,  the  computational 
procedure  corrects  for  boundary-layer  displacement  effects.  With 
the  exception  of  the  boundary-layer  computation,  all  phases  of  the 


procedure  are  accomplished  on  a  TRS-80  Model  II  (z80-based)  micro¬ 
computer.  While  the  boundary-layer  computations  could  also  be 
done  on  the  TRS-80,  run  times  are  too  long  (30  hours)  to  render 
such  runs  useful  in  a  design  environment . 

With  the  exception  of  a  plotting  routine  (SPLOT)  and  an  inter¬ 
polation  routine  (INTERP),  all  programs  have  been  developed  in 
FORTRAN  (ANSI-66)  and,  with  trivial  modifications,  can  thus  be  run 
on  virtually  all  scientific  computers. 

In  Section  2  we  first  present  an  overview  of  the  computational  pro 
cedure.  Next  we  describe  the  various  programs  including  input  and 
output  (program  listings  are  located  in  the  Appendix).  Then,  we 
present  results  of  two  sample  computations. 

Section  3  summarizes  and  discusses  results  of  the  project. 


2 .  ANALYSIS 


2.1  OVERVIEW  OP  THE  COMPUTATIONAL  PROCEDURE 

In  devising  the  computational  procedure,  our  plan  has  been  to 
create  a  method  which  lays  the  foundation  for  a  ship  designer  to 
effect  design  modifications  in  an  interactive  mode.  Because  we 
are  creating  the  foundation,  the  procedure  has  been  developed  in 
several  discrete  modules  each  of  which  performs  a  straightforward 
task. 

Our  starting  point  is  a  disk  file  containing  the  three-dimensional 
hull  coordinates  from  bow  to  stern.  This  file  presumably  has  been 
created  by  the  designer.  As  illustrated  in  Figure  1,  we  first  use 
program  PRMESH  to  generate  a  preliminary  orthogonal,  rectangular 
mesh  for  use  in  computing  the  inviscid  velocities  on  the  hull. 
Before  proceeding  to  computation  of  the  hull  velocities,  program 
SPLOT  can  be  used  to  visually  inspect  the  coordinates  via  the  micro 
computer  video  display.  This  step  provides  an  opportunity  to  spot 
obvious  input  data  errors.  Once  we  are  confident  the  input  data 
are  satisfactorily  smooth,  we  use  program  SHIP  to  compute  the  In¬ 
viscid  velocities  on  the  ship-hull  surface.  This  completes  the 
inviscid  phase  of  the  computational  procedure. 

The  second  (viscous)  phase  of  the  procedure  begins  with  program 
SHPMSH  which  generates  a  nonorthogonal  finite-difference  grid  for 
use  in  the  3-D  boundary-layer  computation,  including  all  required 
geodesic  curvatures  and  metrics.  Then,  using  program  VELOC  we 
interpolate  the  inviscid  velocities  onto  the  nonorthogoral  grid. 

At  this  point,  all  required  input  data  for  program  EDDY3  have  been 
prepared.  The  final  step  in  the  viscous  phase  is  execution  of 
program  EDDY3  which  computes  boundary-layer  development  over  the 
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Figure  1.  Schematic  representation  of  the  computational  procedure 


Upon  completion  of  the  viscous  phase,  the  designer  can  opt  to  in¬ 
clude  effects  of  boundary-layer  displacement.  Should  this  option 
be  selected,  a  disk  file  is  available  from  EDDY3  which  containes 
the  computed  displacement  thicknesses  on  the  nonorthogonal  grid. 
Program  INTERP  interpolates  the  contents  of  this  file  onto  the 
rectangular  grid  used  by  program  SHIP.  The  complete  computational 
sequence  is  then  repeated. 

Certainly  visual  inspection  of  the  inviscid  velocities,  nonorthogon¬ 
al  grid  coordinates,  final  output  from  EDDY3,  etc.  would  be  help¬ 
ful,  if  not  essential,  in  a  design  application.  Simple  modifica¬ 
tions  to  SPLOT,  for  example,  would  provide  the  vehicle  for  such 
inspection.  Even  more  powerful  microcomputer  graphics  techniques 
are  available  for  both  video  and  hard-copy  displays  and  there  seems 
little  point  to  going  to  great  lengths  at  this  point  to  rival  exist¬ 
ing  software.  Thus,  we  have  not  devised  the  needed  modifications 
because  of  the  preliminary  nature  of  this  project  which  has  been, 
to  some  extent,  a  feasibility  study.  Nevertheless,  the  capability 
for  inspection  of  the  input  data  has  been  included  to  give  an  in¬ 
dication  of  the  practicability  of  a  microcomputer-based  approach. 

2.2  PROGRAM  INPUT/OUTPUT 

We  turn  now  to  the  input  and  output  of  the  seven  program  modules. 

The  Appendix  contains  complete  listings  of  each  program.  In  the 
interest  of  clarity,  for  each  program,  we  reiterate  its  overall 
purpose  below. 

2.2.1  Program  PRMESH 

The  purpose  of  this  program  is  to  generate  a  rectangular  finite- 
difference  grid  for  use  in  computing  the  inviscid  velocity  com¬ 
ponents. 

Program  input  data  consists  of  two  disk  files,  viz,  MESH  and  PRMESHIN 
The  first  of  these  two  files  contains  the  coordinates  x,  y,  z  of 
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the  ship  (Figure  2),  where  the  coordinates  are  most  conveniently 
supplied  in  nondimensional  form  (normalized,  for  example,  by  hull 
length,  L).  Note  that  flow  past  the  hul'1  is  toward  positive  x, 
positive  y  is  to  starboard,  and  z  increases  as  we  move  from  keel 
to  waterline. 

In  disk  file  MESH,  the  coordinates  are  presented  on  constant-x 
sections  with  each  line  giving  two  pairs  of  coordinates.  The  for¬ 
mat  is  (3F10. 5,1X,3F10. 5) •  Note  that  the  space  between  coordinate 
pairs  (the  IX)  is  included  for  consistancy  with  the  input  data 
format  specified  for  the  1980  ITTC/SSPA  Workshop  on  Ship  Bouno  y 
Layers.  Thus,  a  few  typical  lines  of  data  assume  the  followir 
form: 

SAMPLE  ’MESH”  INPUT  DATA 


-0. 40000 
-0.40000 
-0. 40000 
-0.35000 
-0.35000 
-O. 35000 


-0. 13152 
-0. 13603 
-0. 13705 
0.00000 
-0.01498 
-0.03952 


-0.07847 
-0.05375 
—0. 01872 
-0.11808 
-0.11790 
-0.11732 


-0.40000 
-0.40000 
—0. 40000 
-0.35000 
-0.35000 
-0.35000 


-O. 13425 
-0. 13668 
—0. 13738 
-0.00690 
-0.02537 
-0.05698 


-0.06825 

-0.03680 

f\  ortftrtft 

•  V  V  V  *W*  V 

-0.11803 

-0.11767 

-0.11668 


Note  that  points  on  each  section  are  input  from  keel  (where  y  =  0) 
to  waterline  (where  z”  -  0).  In  its  current  form,  program  PRMESH 
will  accept  data  on  up  to  50  sections  with  a  maximum  of  20  (y,z) 
pairs  per  section.  One  restriction  on  the  input  data  is  that  an 
even  number  of  (y,z)  pairs  must  be  prescribed.  Note  also  that  the 
input  coordinates  must  define  the  ship  hull  from  bow  to  stern;  a 
partial  mesh  is  insufficient. 

The  second  input  disk  file,  PRMESHIN,  contains  the  following  data: 

IMAX  ■  Number  of  sections  (constant  x)  in  the  computed 
rectangular  mesh  (cannot  exceed  41). 

JMAX  »  Number  of  points  per  section  in  the  computed 
rectangular  mesh  (cannot  exceed  20). 

NXTL  *  Number  of  sections  in  the  input  mesh  of  disk  file 
MESH  (cannot  exceed  50). 

NZTL  *  Number  of  (y,z)  pairs  per  section  in  the  input  mesh 
of  disk  file  MESH  (cannot  exceed  20). 
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X 


a  Specified  streamwise  (x-coordinate)  array  for  the 
computed  rectangular  mesh  (a  total  of  IMAX  values). 

Z  *  Specified  vertical  (z-coordinate)  array  for  the 

computed  rectangular  mesh  (a  total  of  JMAX  values). 

The  input  format  for  this  disk  file  has  been  designed  to  resemble 
NAMELIST  input.  For  IMAX,  JMAX,  NXTL  and  NZTL,  the  format  is 
(1X,3A4,I4)  while  for  X  and  Z  it  is  (1X,3A4 ,E13 • 6) .  A  typical  in¬ 
put  file  is  shown  below: 

SAMPLE  * PRMESHIN  *  INPUT  DATA 


IMAX 

=s 

5 

JMAX 

35 

4 

NXTL 

sr 

40 

NZTL 

ss 

20 

X(l) 

ss 

-0. 500000E 

00 

X  (2) 

s 

-0. 2S0000E 

00 

X  (3) 

=s 

0. OOOOOOE 

00 

X  (4) 

S3 

0. 250000E 

00 

X  (5) 

= 

0. 500000E 

00 

ZU> 

= 

0. OOOOOOE 

CO 

Z  (2) 

* 

-0. 030000E 

00 

Z  (3) 

= 

-0.0&5000E 

00 

Z  (4) 

SS 

-0. 100000E 

00 

As  illustrated,  the  (1X3A4)  part  of  the  format  permits  input  of 
the  variable  name.  Unlike  NAMELIST  input,  all  variables  must  be 
input  and  the  order  of  input  cannot  be  mixed.  An  exception  to  this 
is  that  only  IMAX  (JMAX)  values  of  X(Z)  must  be  input;  the  unused 
parts  of  these  arrays  needn't  be  defined. 

Program  output  is  (1)  a  printed  listing  of  input  disk  file  PRMESHIN 
and  (2)  a  disk  file  named  MESHPR  containing  the  coordinates  of  the 
computed  rectangular  mesh.  Inspection  of  the  latter  can  be  accom¬ 
plished  with  program  SPLOT  described  in  the  next  Subsection. 

2.2.2  Program  SPLOT 

The  purpose  of  program  SPLOT  is  to  permit  section-by-section  visual 
inspection  of  the  rectangular  mesh  generated  by  program  PRMESH. 


Such  an  Inspection  will  aid  in  identifying  smoothness  of  the  mesh 
as  well  as  input  data  errors. 

The  program  is  written  in  BASIC,  makes  use  of  special  features 
of  the  TRS-80  Model  II,  and  cannot  be  readily  adapted  to  another 
machine.  However,  for  reasons  noted  in  the  Introduction,  it  has 
been  included  for  demonstration  purposes  and  is  not  essential  in 
the  overall  computational  procedure.  Input  to  the  program  consists 
of  (1)  prompted  requests  for  IMAX  and  JMAX  and  (2)  disk  file  MESHPR. 
Also,  as  the  visual  inspection  proceeds,  additional  prompted  requests 
are  made  as  will  be  discussed  below. 

All  input  and  output  appears  on  the  video  display.  Upon  execution 
of  SPLOT,  the  program  will  print  "IMAX®?" .  In  response,  the  value 
of  IMAX  used  in  running  PRMESH  must  be  entered.  Next  the  message 
"JMAX=?"  is  displayed.  Again,  the  value  used  in  running  PRMESH 
must  be  entered. 

The  program  will  now  proceed  section  by  section.  On  each  section 

it  will  print  the  streamwise  section  number  and  the  value  of  x  on 

the  section.  This  will  be  followed  by  a  printout  identifying  the 

maximum  and  minimum  values  of  y  and  z  on  the  section.  The  program 

then  requests  the  appropriate  plotting  range,  viz,  the  desired  values 

of  y„_  >  z  .  .  z _ .  After  the  first  section  has  been  in- 

•'mln*  "'max  min  max 

vestigated,  the  program  first  asks  if  the  same  plotting  range  used 
for  the  previous  plot  is  desired.  Entering  "Y”  obviates  the  need 
to  reenter  ym^n»  etc.  Entering  "N”  permits  changing  the  plotting 
range.  This  feature  has  been  included  to  permit  viewing  all  sections 
on  the  same  scale. 

The  plot  for  a  given  section  will  be  displayed  until  a  character 
is  entered.  The  following  options  are  available: 

C  *  Continue  inspection  on  this  section; 

N  *  Next  section  -  proceed  to  same; 

For  any  other  ASCII  character,  the  run  terminates.  The  "C”  option 
permits  changing  the  plotting  range.  This  might  be  done  to  expand 
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the  scale  in  a  particular  region  for  closer  inspection.  Note 
that  if  this  is  done,  the  entered  ym^n»  etc.  become  the  new 
plotting  range  for  all  following  sections.  Hence,  returning  to 
the  original  plotting  range  requires  the  appropriate  input  once 
the  user  proceeds  to  the  next  section.  The  "C"  option  can  be  used 
as  many  times  as  desired  on  each  section. 

As  a  safeguard  against  incorrect  entries,  the  program  asks  if  "ANY 
ERRORS  ?"  have  been  made  each  time  ym^n»  etc.  are  entered.  If  no 
plotting-range  errors  have  been  made,  enter  "N"  and  the  program 
continues.  If  an  error  has  been  made,  enter  "Y"  and  the  correct 
values  can  then  be  typed  in.  Finally,  any  points  which  are  off 
scale  will  be  Indicated  in  the  lefthand  margin  of  the  plot. 

When  the  user  is  satisfied  the  mesh  is  sufficiently  smooth  and  free 
of  errors,  he/she  is  ready  to  proceed  to  computation  of  the  inviscid 
flow.  This  end  is  accomplished  with  the  program  described  in  the 
next  subsection. 

2.2.3  Program  SHIP 

The  purpose  of  program  SHIP  is  to  compute  the  inviscid  velocity 
components  on  the  surface  of  the  ship  hull.  As  an  auxiliary  out¬ 
put,  the  program  computes  the  wave  drag. 

Input  to  program  SHIP  consists  of  two  disk  files,  viz:  (1)  file 
MESHPR  created  by  program  PRMESH;  and  (2)  disk  file  SHIPIN.  The 
latter  file  has  again  been  created  in  a  format  similar  to  NAMELIST. 
The  input  data  are: 

FR  =  Froude  number  defined  by  F  *  U//gL  where  U  is  free- 
stream  velocity,  g  is  gravitational  acceleration  and 
L  is  the  length  scale_used  to  nondimensionalize  the 
ship  coordinates  x,y,z. 

IMAX  »  Number  of  sections  (constant  x)  in  the  rectangular 

mesh  (same  value  used  in  disk  file  PRMESHIN  -  cannot 
exceed  41); 

JMAX  =  Number  of  points  per  section  in  the  rectangular  mesh 
(same  value  used  in  disk  file  PRMESHIN  -  cannot 
exceed  10); 


NINT2  *  Number  of  points  used  in  Chebyshev  quadrature  to 
evaluate  an  integral  whose  dimensionless  wave- 
number  range  is  1  to  •  (cannot  exceed  80); 

MINT3  *  Number  of  points  used  in  Chebyshev  quadrature  to 
evaluate  an  integral  whose  dimensionless  wave- 
number  range  is  zero  to  1  (cannot  exceed  80); 


NZOT  *  Print  flag.  For  zero  the  cell  geometry  and  velocity 
components  are  printed.  Otherwise,  only  the  wave 
drag  is  printed. 

As  with  disk  file  PRMESHIN  of  Subsection  2.2.1,  the  format  is 
(1X,3A4,E13. 6)  for  FR  and  (1X,3A4,I4)  for  all  integer  input  quan¬ 
tities.  A  typical  input  file  is  as  follows 


SAMPLE  ’ SHIPIN ’  INPUT  DATA 


FR 

2 . 500000E-0 1 

I  MAX 

* 

5 

JMAX 

2S 

4 

MINT2 

20 

MINT3 

* 

20 

NZOT 

= 

0 

Note  that  the  maximum  value  of  JMAX  allowed  has  been  limited  to 
10  to  keep  execution  times  on  a  TRS-80  Model  II  within  practical 
limits. 

Output  from  program  SHIP  consists  of:  (1)  video  messages  indicating 
when  each  of  three  stages  of  the  computation  is  completed;  (2)  print 
ed  output  depending  on  the  value  assigned  to  input  flag  NZOT;  and 
(3)  a  disk  file  named  UEWE  containing  mesh  coordinates  and  computed 
velocity  components. 

The  first  part  of  the  printout  is  diskfile  SHIPIN.  This  provides 
an  immediate  check  that  input  data  have  been  correctly  entered. 


Printed  messages  during  the  computation  appear  on  the  video  display 
upon  completion  of  calls  to  subroutines  VEL0C1,  VEL0C2  and  VEL0C3. 

The  message  after  completing  a  call  to  VELOCn  is  "VELOCn  COMPUTATIONS 


COMPLETE" . 
of  the  run 


These  messages  permit  the  user  to  monitor  the  progress 


The  next  part  of  the  printout  depends  upon  the  value  assigned  to 
NZOT.  If  NZOT  is  zero  the  program  prints:  the  three  coordinates 
of  the  rectangular  cell  centroids;  the  computed  dimensionless 
velocity  components  u/U,  v/U,  w/U;  and  the  dimensionless  total 
velocity  q/U.  The  print  proceeds  section  by  section  including 
incremental  wave  drag  on  the  section  and  total  wave  drag  computed 
up  to  the  section. 

The  final  part  of  the  printout  is  the  total  dimensionless  wave 
drag  (RWAVE)  defined  in  terms  of  the  drag  R  by 

RWAVE  =  — ^ -  (1) 

%pU2L2 

where  p  is  density,  while  U  and  L  are  as  defined  above. 

Finally,  dlskflle  UEWE  is  written  for  use  by  program  VELOC  described 
in  Subsection  2.2.5.  Upon  completion  of  this  phase  of  the  computa¬ 
tion,  we  have  determined  the  invlscid  velocity  components  on  an 
orthogonal  rectangular  mesh. 

2.2.4  Program  SHPMSH 

The  purpose  of  this  program  is  to  generate  a  nonorthogonal  finite- 

difference  mesh  for  use  by  the  three-dimensional  boundary-layer 

program  EDDY3  (Subsection  2.2.6).  The  program  generates  least- 
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squares  cubic-spline  fits  to  accurately  and  smoothly  compute  all 
required  geodesic  curvatures  and  metric  coefficients. 

Input  data  for  SHPMSH  consists  of  two  dlskflles,  viz,  MESH  and 
NAMLIST.  Diskfile  MESH  is  the  original  file  containing  the  ship- 
hull  coordinates.  In  the  interest  of  economy,  this  file  needn't 
cover  the  entire  hull.  Rather,  it  need  only  cover  the  part  of  the 
hull  over  which  the  boundary  layer  is  to  be  computed.  Diskfile 


NAMLIST  contains  the  following  input  data: 

ISYM  ■  Symmetry  flag.  To  force  mesh  symmetry  (dy/dz=0) 

at  the  waterline  use  ISYM=0.  Otherwise  the  program 
automatically  infers  the  appropriate  value  of  dy/dz 
at  the  waterline. 

IU1  *  Logical  unit  number  for  input  diskfile  MESH. 

IU2  =  Logical  unit  number  for  temporary  diskfile  SCRATCH1. 

IU3  *  Logical  unit  number  for  temporary  diskfile  SCRATCH2. 

IU4  *  Logical  unit  number  for  output  diskfile  GEOMETRY. 

IU5  *  Logical  unit  number  for  final  printout. 

IMAX  *  Number  of  sections  in  the  input  mesh  of  diskfile  MESH 
(cannot  exceed  40). 

JMAX  =  Number  of  (y,¥)  pairs  per  section  in  the  input  mesh 
of  diskfile  MESH  (cannot  exceed  40). 

NXTL  »  Number  of  sections  in  the  computed  nonorthogonal 
mesh  (cannot  exceed  40). 

NZTL  ■  Number  of  points,  per  section  in  the  nonorthogonal 
mesh  (cannot  exceed  21). 

JFZ  *  Number  of  key  points  used  for  the  least-squares 

cubic-spline  fit  in  the  girthwise  direction  (cannot 
exceed  30). 

LZ  =  Array  of  indices  of  the  key  points  to  be  used  in  the 
girthwise  direction. 

JFX  *  Number  of  key  points  used  for  the  least-square  cubic- 
spline  fit  in  the  streamwise  direction  (cannot 
exceed  30). 

LX  =  Array  of  indices  of  the  key  points  to  be  used  in  the 
streamwise  direction. 

PSIZ  =  Nondimensional  distance  (measured  from  the  bow)  of 

the  most  upstream  section  in  the  nonorthogonal  mesh. 

XSTART  *  Value  of  x  at  the  most  upstream  section  in  the 
nonorthogonal  mesh. 

XEND  =  Value  of  x  at  the  most  downstream  section  in  the 
nonorthogonal  mesh. 

DX  *  Array  of  streamwise  stepsizes,  Ax^,  for  the  nonorthogon 
al  mesh  (cannot  exceed  10). 

NX  *  Array  of  indices  indicating  the  number  of  times  Ax^ 
is  repeated,  viz,  the  mesh  will  have  Ax=DX(l)  for 
the  first  NX(1)  stations,  Ax=DX(2)  for  the  next  NX(2) 
Stations,  etc. 


DZ  ■  Array  of  girthwise  stopsizes,  Az, ,  for  the  non- 

orthogonal  mesh  (cannot  exceed  ID).  Note  that  the 
mesh  is  constructed  so  that  z=0  on  the  keel  and  z=l 
on  the  waterline. 

NZ  ■  Array  of  indices  indicating  the  number  of  times 
A7j  is  repeated. 

As  in  programs  PRMESH  and  SHIP,  the  format  for  integer  quantities 
is  (1X,3A4,I4)  while  floating  point  variables  use  (1X,3A4,E13.6) . 
Sample  input  is  as  follows: 


SAMPLE  'NAMLIST*  INPUT  DATA 


ISYM 

= 

1 

DX  (5) 

8 

O.OOOOOOE 

00 

IU1 

32 

7 

DX  (6) 

8 

O.OOOOOOE 

00 

IU2 

32 

8 

DX  (7) 

8 

O.OOOOOOE 

00 

IU3 

32 

9 

DX  (8) 

8 

O.OOOOOOE 

00 

IU4 

23 

10 

DX  (9) 

8 

O.OOOOOOE 

00 

IU5 

32 

2 

DX ( 10) 

S 

O.OOOOOOE 

00 

I  MAX 

= 

30 

NX  ( 1 ) 

8 

5 

JMAX 

32 

18 

NX  (2) 

8 

2 

NXTL 

32 

13 

NX  (3) 

8 

7 

NZTL 

* 

11 

NX  (4) 

8 

0 

JFZ 

8 

NX  (5) 

32 

0 

LZ  ( 1  > 

32 

1 

NX  (6) 

8 

0 

LZ<2> 

= 

3 

NX  (7) 

8 

0 

LZ  (3) 

32 

5 

NX  (8) 

8 

0 

LZ(4) 

32 

7 

NX  (9) 

8 

0 

LZ  (3) 

S 

10 

NX (10) 

8 

0 

LZ  (6) 

8 

12 

ZSTART 

= 

O.OOOOOOE 

00 

LZ  (7) 

= 

15 

ZEND 

8 

1 .  OOOOOOE 

00 

LZ  (8) 

a 

18 

DZ  ( 1 ) 

8 

1 . OOOOOOE- 

-01 

JFX 

* 

12 

DZ  (2) 

8 

0.  OOOOOOE 

00 

LX  ( 1 ) 

= 

1 

DZ  (3) 

= 

O.OOOOOOE 

00 

LX  (2) 

= 

3 

DZ  (4) 

8 

O.OOOOOOE 

00 

LX  (3) 

* 

6 

DZ  (5) 

32 

O.OOOOOOE 

00 

LX  (4) 

32 

9 

DZ  (6) 

8 

O.OOOOOOE 

00 

LX  (5) 

= 

12 

DZ  (7) 

8 

O.OOOOOOE 

00 

LX  (6) 

3 

15 

DZ  (8) 

8 

O.OOOOOOE 

00 

LX  (7) 

32 

18 

DZ  (9) 

8 

O.OOOOOOE 

00 

LX  (8) 

8 

21 

DZ ( 10) 

8 

0. OOOOOOE 

00 

LX  (9> 

as 

24 

NZ  ( 1 ) 

32 

10 

LX ( 10) 

33 

26 

NZ  (2) 

8 

0 

LX  (11) 

33 

28 

NZ  (3) 

8 

0 

LX (12) 

32 

30 

NZ  (4) 

8 

0 

PSIZ 

8 

3. OOOOOOE-Ol 

NZ  (5) 

8 

0 

XSTART 

8 

-7.000000E-01 

NZ  (6) 

a 

0 

XEND 

8 

9. 000000E-01 

NZ  (7) 

8 

0 

DX(1) 

32 

1 . OOOOOOE-Ol 

NZ  (8) 

8 

0 

DX  (2) 

S 

1 . 500000E-01 

NZ  (9) 

8 

0 

DX  (3) 

8 

1. OOOOOOE-Ol 

NZ ( 10) 

8 

0 

DX  (4) 

8 

O.OOOOOOE  00 

Output  from  program  SHPMSH  consists  of:  video  display  messages 
indicating  when  various  phases  of  the  computation  are  complete; 
printed  output  of  computed  quantities;  and  diskfile  GEOMETRY 
for  use  by  programs  VELOC  (Subsection  2.2.5)  and  EDDY3  (Subsection 
2.2.6). 

The  first  printed  output  is  a  listing  of  diskfile  NAMLIST  which 
permits  checking  for  input  data  errors. 

Three  messages  are  then  printed  on  the  video  display  during  the 
run  to  permit  monitoring  the  progress  of  the  run.  The  messages 
are: 

1.  "All  Section  Data  Smoothed  and  Differentiated"; 

2.  "Geometrical  Parameters  Computed  for  J=j"; 

3*  "d(Theta)/dz  Computed  for  All  Sections". 

The  second  message  is  repeated  for  each  girthwise  station  where  j 
varies  from  1  to  JMAX. 


The  next  output  is  the  computed  quantities  along  each  girthwise 
station.  This  information  is  routed  to  the  printer  if  IU5=2  or 
saved  on  disk  provided  IU5  is  greater  than  or  equal  to  7 •  The 
data  is  saved  on  diskfile  F0RTu5/DAT  where  "u5"  is  the  numerical 
value  of  IU5  when  IU53.7-  This  option  permits  reiterating  to  in¬ 
clude  displacement  effects  via  program  INTERP  (Subsection  2.2.7) 

The  final  output  is  diskfile  GEOMETRY  which  is  used  by  programs 
VELOC  and  EDDY 3- 


2.2.5  Program  VELOC 

The  purpose  of  this  program  is  to  interpolate  the  inviscid  velo¬ 
cities  onto  the  nonorthogonal  coordinate  mesh  for  use  in  the  three 
dimensional  boundary-layer  computation. 

Input  data  for  this  program  consists  of  three  diskfiles,  viz, 
NAMVEL,  UEWE  and  GEOMETRY.  Diskfiles  UEWE  and  GEOMETRY  are  creat¬ 
ed  by  programs  SHIP  and  SHPMSH,  respectively.  Diskfile  NAMVEL 


contains  the  following  input  data: 


ISYM  *  Symmetry  flag.  To  treat  the  waterline  as  a  symmetry 
plane,  use  ISYM=0.  Otherwise,  the  program  will 
automatically  infer  velocities  at  the  free  surface. 

IU1  =  Logical  unit  number  for  input  diskfile  UEWE. 

IU2  *  Not  used. 


IU3  =  Logical  unit  number  for  output  diskfile  VELOCITY. 

IU4  =  Logical  unit  number  for  input  diskfile  GEOMETRY. 

IU5  =  Logical  unit  number  for  final  printout . 

IMAX  =  Number  of  sections  in  the  input  mesh  of  diskfile 
UEWE  (cannot  exceed  40). 

JMAX  =  Number  of  points  per  section  in  the  input  mesh  of 
diskfile  UEWE  (cannot  exceed  19)* 

NXTL  =  Number  of  sections  in  the  nonorthogonal  mesh  of 
diskfile  GEOMETRY  (cannot  exceed  40). 

NZTL  =  Number  of  points  per  section  in  the  nonorthogonal 
mesh  of  diskfile  GEOMETRY  (cannot  exceed  21). 

JFZ  *  Number  of  key  points  for  the  least-squares  cubic- 

spline  fit  in  the  girthwise  direction  (cannot  exceed  30). 

LZ  *  Array  of  indices  of  the  key  points  to  be  used  in 
the  girthwise  direction. 

JFX  =  Number  of  key  points  used  for  the  least-squares  cubic- 
spline  fit  in  the  streamwise  direction  (cannot 
exceed  30). 

LX  *  Array  of  indices  of  the  key  points  to  be  used  in 
the  streamwise  direction. 

Note  that  because  program  SHIP  does  not  provide  keel  or  waterline 
velocities,  this  program  computes  both.  In  so  doing,  the  program 
replaces  JMAX  by  JMAX+2.  Also,  JFZ,  JFX,  LZ  and  LX  needn't  be 
the  same  as  used  in  program  SHPMSH.  As  in  all  similar  input  data 
files  of  preceding  programs,  the  format  is  (1X,3A4,I4)  for  integer 
input  quantities  and  (1X,3A4,E13. 6)  for  floating  point  data.  A 
typical  listing  of  diskfile  NAMVEL  follows. 
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SAMPLE  ’NAMVEL’  INPUT  DATA 


ISYM 

s 

1 

LZ  (6) 

SS 

13 

IU1 

ss 

7 

LZ  (7) 

= 

16 

IU2 

ss 

8 

LZ  (8) 

= 

19 

IU3 

= 

9 

JFX 

= 

12 

IU4 

ss 

10 

LX  ( 1 ) 

ss 

1 

IU5 

= 

2 

LX  (2) 

= 

3 

WAX 

ss 

31 

LX  (3) 

= 

6 

JMAX 

= 

17 

LX  (4) 

= 

9 

NXTL 

ss 

15 

LX  (5) 

ss 

12 

NZTL 

ss 

11 

LX  (6) 

s 

15 

JFZ 

ss 

8 

LX  <7) 

ST 

18 

LZ  ( 1 ) 

- 

1 

LX  (8) 

= 

21 

LZ  (2) 

ss 

3 

LX  <9) 

= 

24 

LZ  (3) 

ss 

5 

LX (10) 

as 

27 

LZ  (4) 

=: 

7 

LX  (11) 

= 

29 

LZ  (5) 

St 

10 

LX (12> 

ss 

31 

Program  output  consists  of:  video  messages  indicating  when 
specific  phases  of  the  computation  are  complete;  printed  output 
listing  interpolated  nonorthogonal  velocity  components;  and 
diskfile  VELOCITY  for  use  by  program  EDDY3. 

The  first  output  is  a  printout  of  input  diskfile  NAMVEL  for  in¬ 
spection  to  insure  no  input  data  errors  have  been  made. 

The  next  output  phase  consists  of  the  following  messages  on  the 
video  display: 

1.  ’’Section  Data  Smoothed  &  Differentiated  for  I=i"; 

2.  ’’All  Section  Data  Smoothed  &  Differentiated”. 

The  first  message  is  displayed  for  each  section  where  i  ranges 
from  1  to  NXTL. 
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The  program  then  provides  printed  output  for  all  computed 
velocity  components  and  for  dw/dz  on  the  keel  and  waterline. 

Finally,  dlskflle  VELOCITY  Is  written  for  use  by  program  EDDY3* 


2.2.6  Program  EDDY3 

The  purpose  of  this  program  is  to  compute  three-dimensional 
boundary-layer  development  on  the  ship  hull.  While  all  of  the 
other  programs  listed  in  the  Appendix  are  shown  in  their  forms 
appropriate  to  the  TRS-80  Model  II  microcomputer,  the  listing  of 
EDDY3  is  appropriate  to  a  UNIVAC  11 08  computer.  Required  changes 
for  CDC  or  IBM  computers  are  confined  to  SUBROUTINE  STORIT  and  are 
Included  via  Comment  lines. 

Input  data  for  program  EDDY3  consists  of:  diskfile  GEOMETRY; 
diskfile  VELOCITY;  three  NAMELIST  files,  viz,  NAME,  DATA  and 
STRT;  and  two  lines  of  card-type  data.  Diskfiles  GEOMETRY  and 
VELOCITY  are  created  by  programs  3HPMSH  and  VELOC,  respectively. 

The  rest  of  the  input  is  described  below. 

The  first  line  of  input  is  read  from  unit  TAPEIN  which  has  a 
default  value  of  5-  This  input  is  the  title  of  the  case  to  be 
run  which  can  be  up  to  5^  characters  long  (format  9A6). 

The  next  segment  of  the  input  is  also  read  from  unit  TAPEIN  and 
is  in  NAMELIST  format.  The  data  are  as  follows: 

NAMELIST/NAME/ :  This  namelist  includes  primary  parameters  defining 

the  grid,  print  logic,  etc. 


DETA  3  Array  of  An  values;  only  the  An  nearest  the  surface 
need  be  input.  Note  that  n  is  the  conventional 
Levy-Lees-type  normal  coordinate. 

EPS  *  Convergence  criterion  for  turbulence-model  equations 
(DEFAULT3 . 01) . 

EPST  3  Convergence  criterion  for  crossflow  momentum  equation 
(DEFAULT3. 01). 

EPSV  3  Convergence  criterion  for  streamwise  momentum  equation 
(DEFAULT3. 01). 


ETAE  = 

ICHORD  * 
IPLOW  = 


IPPRNT  * 


IPX  * 

IPZ  * 

ISPAN  = 
ITMAX  = 
NPT  = 

ATfPO  = 

i* 

NXSTRT  = 
NXT  = 
NZSTRT  * 
NZT  * 

TAPEIN  = 

TAPEOT  * 
TAPEGP  = 

TAPEPF  * 

TAPEDT  - 

TAPEVL  = 

VGP  « 
ISHORT  = 


Boundary-layer  edge  value  for  n  at  Initial  station 
(DEFAULTS). 

Not  used  ( DEFAULT* 1) . 

Flag  indicating  crossflow  integration  direction. 
IFL0W=1  for  integration  from  keel  to  waterline  and 
IFL0W=2  for  integration  from  waterline  to  keel 
( DEFAULT* 1) . 

Print  flag.  IFPRNT=0  deletes  printout  of  contents 
of  diskfiles  GEOMETRY  and  VELOCITY.  Otherwise  the 
contents  of  these  two  input  files  are  printed 
( DEFAULT* 0) . 

Intervals  at  which  profiles  are  printed  in  stream- 
wise  direction,  e.g.,  IPX=10  means  print  every 
tenth  step. 

Intervals  at  which  profiles  are  printed  in  girthwise 
direction. 

Not  used  ( DEFAULT* 1) . 

Maximum  number  of  iterations  allowed  (DEFAULT=20) . 

Maximum  number  of  mesh  points  allowed  normal  to 
hull  (must  not  exceed  101...  DEFAULT=101) . 

Lasiinar/turbulent  flow  fists.  To  st2.i*tr  f r* cm  1  stminsti* 
profiles  use  NTR=0.  Otherwise,  start  from  turbulent 
profiles  (DEFAULT=1) . 

Index  of  initial  section,  (must  be  1). 

Index  of  final  section  to  which  computation  proceeds. 

Index  of  initial  girthwise  line  (Must  be  1). 

Index  of  final  girthwise  line  to  which  computation 
proceeds. 

Logical  unit  number  for  subsequent  NAMELIST  input 
( DEFAULT*  5 ) • 

Logical  unit  number  for  printed  output  (DEFAULT=6). 

Logical  unit  number  for  input  diskfile  GEOMETRY 
(DEFAULT*l6) . 

Logical  unit  number  for  temporary  diskf ile-this  is 
a  very  large  file  (DEFAULT=17) . 

Logical  unit  number  for  temporary  diskfile 
(DEFAULT*l8) . 

Logical  unit  number  for  input  diskfile  VELOCITY 
(DEFAULT*19). 

Geometric  progression  ratio  for  normal  grid. 

Flag  for  slightly  abbreviated  print.  ISH0RT=0 
deletes  details  of  convergence  parameters.  Other¬ 
wise  all  details  are  given  (DEFAULT*0). 
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RFTRB 


Relaxation  coefficient  for  turbulence  model  equations 
(DEFAULT®. 8) . 

RFVEL  =  Relaxation  coefficient  for  momentum  equations 
( DEFAULT® 1) . 

NAMELIST/DATA/ :  This  namelist  defines  the  various  turbulence- 
model  closure  coefficients. 


SALFA  = 
SALFAS  = 
SBETA  = 
SBETAS  = 
SLAMDA  = 
SSIGMA  = 
SSGMAS  = 
USTOP  = 

ZIOTAE  = 

ZIOTAL  = 

RW2  * 


The  closure  coefficient  a  (DEFAULT-1.11111). 

The  closure  coefficient  a*  (DEFAULT® . 3) . 

The  closure  coefficient  3  (DEFAULT® .15) . 

The  closure  coefficient  3*  (DEFAULT® . 09) . 

The  closure  coefficient  X  (DEFAULT® . 091) . 

The  closure  coefficient  a  ( DEFAULT® . 5 ) . 

The  closure  coefficient  a*  (DEFAULT®. 5) 

The  value  of  u+=u/u  below  which  turbulent  dissipation 
rate  is  analytically  prescribed  (DEFAULT®**). 

Dimensionless  value  of  turbulent  mixing  energy  at 
the  boundary-layer  edge  (DEFAULT=3 . 75 • 10  5). 

Ratio  of  turbulent  length  scale  to  boundary  layer 
thickness  (DEFAULT®. 09) . 

The  closure  coefficeint  R^  (DEFAULT® 4 . ) . 


NAMELIST/STRT/ :  This  namelist  provides  input  parameters  needed 


to  define  turbulent  starting  profiles. 


ALAMI 

CFXI 
CFZI 
DELTA I 
THETXI 
THETZI 


Array  of  values  at  each  girthwise  station  for  ratio 
of  turbulent  mixing  energy  to  the  corresponding 
flat-plate  value. 

Array  of  streamwise  skin  friction,  Cf  ®tw  /(hpU2), 
at  each  girthwise  station.  x 

Array  of  crossflow  skin  friction,  Cf  =tw  /(*spU2) , 
at  each  girthwise  station. 

Array  of  boundary-layer  thickness  at  each  girth¬ 
wise  station. 

Array  of  streamwise  momentum  thickness  at  each 
girthwise  station. 

Array  of  crossflow  momentum  thickness  at  each  girth¬ 
wise  station  (NOTE:  This  array  must  be  the  same  as 
THETXI  in  the  present  form  of  the  program.). 


The  final  Input  Is  a  line  defining  the  freestream  velocity,  static 
pressure  and  kinematic  viscosity.  These  quantities  are  dimensional 
and  the  format  is  (3E12.4). 

Printed  program  output  is  mostly  self  explanatory  and  will  thus 
be  described  only  briefly  here.  The  first  part  of  the  printout 
lists  all  input  parameters.  Then,  on  every  section,  key  integral 
parameters  such  as  cfx,  Cfz,  <5X*,  <5Z*,  0X>  0z  are  6iven  at  all 
girthwise  stations.  At  the  specified  IPX  and  IPZ  intervals,  velocity 
and  turbulence  parameter  profiles  are  printed. 

The  program  also  writes  a  special  diskfile  (use  logical  unit 
number  9)  containing  computed  displacement  thickness.  This  disk- 
file  is  used  in  program  INTERP  (next  Subsection)  to  allow  re¬ 
peating  the  computation  with  displacement  effects  included. 

2.2.7  Program  INTERP 

The  purpose  of  this  program  is  to  interpolate  the  computed  dis¬ 
placement  thickness  distribution  onto  the  orthogonal  rectangular 
coordinate  system  used  by  program  SHIP.  The  whole  computational 
procedure  can  then  be  repeated  to  account  for  displacement  effects. 

The  program  is  written  in  BASIC  and  uses  an  extremely  efficient 
MACHINE  LANGUAGE  interpolation  routine  referred  to  as  USR1  in 
the  program. 

The  first  two  Inputs  to  the  program  are  prompted  by  the  video 
messages  "NXTL=?"  and  "NZTL=?";  the  values  of  NXTL  and  NZTL  used 
in  running  EDDY3  must  be  entered.  The  next  input  is  prompted  by 
the  message  "FILE  NAME  FOR  NONORTHOGONAL  MESH  ?";  the  response  is 
the  name  given  to  the  print  file  generated  by  program  SHPMSH, 
viz,  F0RTu5/DAT.  The  option  is  included  to  enter  the  file  name 
in  case  the  user  opts  to  rename  the  file  after  running  program 
SHPMSH.  The  final  input  is  prompted  by  the  messages  "IMAX*?"  and 
"JMAX*?";  the  values  of  IMAX  and  JMAX  used  in  running  program 
PRMESH  must  be  used. 
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Program  output  is  a  revised  diskfile  MESHPR  containing  the  modifi¬ 
ed  ship-hull  coordinates.  The  original  MESHPR  is  renamed  MESHPR/ 

OLD.  The  contents  of  the  new  MESHPR  diskfile  can  be  examined 
using  program  SPLOT  (Subsection  2.2.2).  Then,  the  entire  computation¬ 
al  sequence  is  repeated  to  determine  the  effect  of  including  dis¬ 
placement  thickness  in  the  inviscid-flow  computation. 

2.3  APPLICATIONS 

In  this  section  we  exercise  the  computational  package  to  exhibit 
accuracy  and  capabilities  for  two  test  cases.  We  first  compute 
the  inviscid  flow  past  a  very  simple  hull  investigated  by  Michell^-. 
Then,  we  do  both  inviscid  and  viscous  computations  for  the  SSPA  720 
ship  hull3. 

2.3-1  Michell* s  Example 

Our  first  test  case  is  the  sinusoidal  hull  form  investigated  by 
Michell.  The  hull  coordinates  are  defined  by 

y  =  -  .02  [l  +  cos  (2ttx)][i  +  cos  (IOttz)]  (2) 

Michell  found  that  for  a  Froude  number  based  on  hull  length  of 
0.25,  the  wave-drag  coefficient  is 

Cp  =  - - —  =  0.60-10-4  (3) 

R  *spU2L2 

In  our  computations  we  use  10  equally-spaced  mesh  points  in  the 
vertical  direction.  The  number  of  points  used  in  the  streamwise 
direction  has  been  varied  to  determine  the  number  needed  to  achieve 
acceptable  engineering  accuracy.  Again  using  equally-spaced 
points,  computations  have  been  made  for  11,  21  and  ^1  points. 

Table  1  summarizes  total  wave-drag  coefficient. 


No.  of  Streamwise 

10,,CR 

Mesh  Points 

11 

0.49 

21 

0.52 

41 

0.56 

Michell 

0.60 

Figure  3  provides  an  interesting  view  of  the  rate  of  convergence 
of  the  method.  The  figure  shows  the  local  drag  coefficient  as  a 
function  of  distance  along  the  hull.  Note  that  the  peak  values 
near  2x/L  *  -0.25  and  2x/L  =  +0.25  are  almost  an  order-of-magnitude 
larger  than  C^. 

As  a  test  of  the  sensitivity  to  vertical-mesh-point  spacing,  we 
have  repeated  the  11-streamwise-point  computation  with  only  6  points 
in  the  vertical  direction.  The  predicted  is  1.3#  lower  than 
the  value  listed  in  Table  1. 

Based  on  results  obtained  in  these  computations,  accuracy  accept¬ 
able  for  preliminary  design  (10#)  can  be  obtained  on  a  mesh  of 
about  15  streamwise  points  by  6  vertical  points.  Such  a  run  takes 
about  45  minutes  on  a  TRS-80  Model  II  microcomputer. 

2.3.2  SSPA  Model  720 

Turning  now  to  the  SSPA  Model  720,  we  exercise  both  the  inviscid 
and  the  viscous  segments  of  the  computational  package.  However, 
for  reasons  which  become  evident  below,  we  do  not  exercise  the  dis¬ 
placement-thickness  correction  phase  for  this  hull. 

Figures  4  and  5  show  computed  streamwise  velocity  distribution  near 
the  free  surface  and  along  the  keel,  respectively.  In  all  cases 
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Figure  3*  Variation  of  predicted  wave-drag  distribution  with  streamwise  mesh 
point  number,  N  . 


Figure  5.  Computed  velocity  distribution 


the  finite-difference  grid  consists  of  41  points  in  the  stream- 
wise  direction  and  11  points  girthwise.  Each  figure  shows  velocity 
distributions  for  Froude  numbers  based  on  total  hull  length  of 
0.15,  0.30  and  0.60.  Also  shown  in  the  figures  are  corresponding 
velocity  distributions  obtained  from  a  Douglas-Neumann  "double- 
hull”  computation. 

Near  the  waterline  the  Fr  =  0.15  distribution  is  quite  oscillatory, 
a  feature  characteristic  of  thin-ship  theory  in  the  limit  of  small 
Froude  number.  In  the  region  between  2x/L  =  -0.30  and  2x/L  =  0.20, 
the  Fr  =  0.30  distribution  has  a  strong  adverse  gradient  which 
almost  certainly  would  lead  to  boundary-layer  separation.  For  Fr  * 
0.60,  the  near-waterline  distribution  corresponds  to  flow  accelera¬ 
tion  over  the  first  70%  of  the  hull.  In  contrast  to  all  of  the 
thin — ship-theory  predictions,  the  "double-hull"  computation  in¬ 
dicates  a  more-or-less  symmetrical  distribution  with  nearly  constant 
velocity  over  the  central  40%  of  the  hull.  Figure  5  indicates 
similar  trends  along  the  keel  with  the  exception  that  no  oscillations 
appear  for  the  Fr  =  0.15  computation. 

Because  the  Fr  *  0.60  distributions  appear  to  be  the  least  likely 
to  cause  numerical  difficulty  in  running  EDDY3,  we  have  selected 
Fr  *  0.60  as  most  appropriate  for  the  test  case  (the  actual  value 
is  not  available  in  our  literature  on  this  hull). 

Figure  6  compares  computed  and  measured  momentum  thickness,  0. 

The  figure  also  shows  computed  results  obtained  using  Douglas- 
Neumann  velocities.  Lines  A,  B  and  C  lie  along  the  keel,  about 
15%  of  the  way  between  keel  and  waterline,  and  about  midway  between 
keel  and  waterline,  respectively.  As  shown,  the  most  notable 
difference  is  the  rapid  increase  in  9  along  Lines  A  and  B  downstream 
of  2x/L  *  0.20.  We  predict  boundary-layer  separation  at  2x/L  =  0.36 
over  a  significant  portion  of  the  hull.  Although  this  result  at 
first  appears  inconsistent  with  the  velocity  distributions  of  Figures 
4  and  5,  there  is  no  contradiction.  The  predicted  crossflow  veloc¬ 
ities  are  quite  large  (w/U~0.25)  and  along  streamlines  we  indeed 
predict  relatively  large  adverse  gradients  approaching  2x/L  *  0.36. 
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Because  the  numerical  predictions  fail  to  cover  such  a  large  region 
of  the  hull,  there  is  little  point  in  reiterating  to  account  for 
displacement  effects  at  this  time. 


3.  DISCUSSION 


The  foundation  has  been  developed  for  a  microcomputer-based  inter¬ 
active  design  procedure  for  general  ship  hulls.  Typical  computing 
time  for  all  but  the  boundary-layer  computation  are  of  the  order 
of  3  to  4  hours  on  a  TRS-80  Model  II  microcomputer.  Using  a 
16-bit  microprocesser,  these  computing  times  would  be  reduced  to 
less  than  a  half  hour.  Typically,  EDDY3  runs  in  about  20  minutes 
on  a  UNIVAC  1108  which  translates  to  about  3  hours  on  a  16-bit 
microcomputer.  Thus,  the  computational  package  constitutes  a 
practical  "numerical  towing  tank"  which  can  be  used  to  effect  and 
test  design  modifications  in  a  single  workday. 

On  the  one  hand,  the  weakest  link  in  the  computational  package  is 
program  SHIP  which  determines  the  inviscid  velocities.  Based  on 
classical  thin-ship  theory.  Its  limitations  are  well  known.  On 
the  other  hand,  program  EDDY3  embodies  one  of  the  most  tested  and 
validated  turbulence  models  available.  The  model  has  even  been 
shown  to  accurately  predict  properties  of  "thick"  boundary  layers 
typical  of  those  near  the  stern  of  a  ship  hull^. 

Research  in  the  immediate  future  should  focus  upon  improving  the 
ability  to  predict  the  inviscid  velocities.  The  formulation  of 

7 

Noblesse  and  Dagan  ,  which  amount  to  a  coordinate  transformation 
using  the  thin-ship  velocity  potential,  should  be  incorporated  in 
program  SHIP.  Also,  normal  pressure  gradient  in  the  boundary-layer 
computation  should  be  accounted  for  to  provide  accurate  predictions 
when  the  boundary  layer  becomes  "thick".  These  two  steps  alone 
would  yield  a  computational  package  which  could  be  immediately 
applied  by  serious  hull  designers.  To  further  aid  the  designer,  a 
master  program  should  be  devised  to  coordinate  input  data  prepara¬ 
tion  for  the  various  program  modules.  Additionally,  available  graphic 
packages  should  be  adapted  to  the  computational  procedure  to  per¬ 
mit  more-detailed  display  and  examination  of  calculated  mesh  and 
flow  details. 


Long  range  research  topics  should  include:  (a)  numerical  studies 
to  aid  development  of  a  more-accurate  integral-method  treatment 
of  the  boundary  layers,  and  (b)  development  of  an  approximate 
engineering-oriented  procedure  for  treating  hulls  on  which  the 
boundary  layer  separates.  The  modular  structure  of  the  computation 
al  package  developed  in  this  project  is  readily  adaptable  to  such 
advanced  developments. 
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APPENDIX:  PROGRAM  LISTINGS 


This  Appendix  contains  listings  of  the  seven  program  modules  which 
constitute  the  computational  package.  In  the  listings,  the  number 
zero  appears  with  a  slash  to  distinguish  it  from  the  letter  ”0” . 
The  various  programs  are  located  on  the  following  pages: 


PROGRAM  NAME  PAGES 

PRMESH . 34-35 

SPLOT . 36-38 

SHIP . 39-45 

SHPMSH . 46-55 

VELOC . . . 56-62 

EDDY  3 . 63-96 

iwTEnr . . . .  . . .  .  .  -  •  yt 


00100 

00200 

00300 

00400 

00300 

00600 

00700 

00000 

00900 

01000 

01100 

01200 

01300 

01400 

01300 

01600 

01700 

01800 

01900 

02000 

02100 

02200 

02300 

02400 

02500 

02600 

02700 

02800 

02900 

03000 

03100 

03200 

03300 

03400 

03500 

03600 

03700 

03800 

03900 

04000 

04100 

04200 

04300 

04400 

04500 

04600 

04700 

04800 

04900 

05000 

03100 

05200 

05300 

03400 

03500 

05600 

05700 

03800 

03900 

06000 


PROGRAM  PRMESH 

C  *************************************************************** 

C  CREATE  RECTANGULAR  MESH  FOR  USE  IN  PROGRAM  SHIP 
C  *************************************************************** 
DIMENSION  XB<50> i Y< 41*20) » YB<50»20>  * YI (50*20) i ZB (20) 
COMMON/DATA/  X<41 ) »  Z(20) 

COMMON/NDATA/  I MAX , JMAXf  NXTL »  NZTL 
CALL  NAM IN 
NZTLM=NZTL-1 
DO  10  1=1 t NXTL 

READ <9* 100)  XB < I ) t  YB (It  NZTL ) t  ZB (NZTL) t YB ( 1 t NZTLM ) t ZB(NZTLM) 

DO  1  J=3t  NZTLM i 2 

JJ-NZTL-J+1 

JM«J-1 

1  READ(9t 110)  YB( I f JJ) , ZB( JJ) * YB( I* JM) *  ZB( JM) 

JN»1 

JS«2 
ZN»ZB( 1 ) 

ZS=ZB(2> 

DO  5  J=1»JMAX 

2  IF(Z(J).LT.ZS.OR.ZN.EQ.ZS)  GO  TO  3 

YI < I t J>=YB( 1 1 JN>  +  ( YB< I t  JS ) —YB ( 1 1  JN> ) * ( Z ( J ) -ZN ) / ( ZS-ZN) 

GO  TO  5 

3  JN**JN+ 1 

IF( JN.GE.NZTL)  GO  TO  4 
JS* JS+ 1 
ZN=ZB( JN> 

ZS=*ZB  <  JS ) 

GO  TO  2 

4  YI(ItJ)=0. 

5  CONTINUE 

10  CONTINUE 
I  W»1 
IE*2 

XW=XB( 1 ) 

XE»XB(2) 

DO  15  1  =  1 1 I MAX 

11  IF(X(I).GT.XE.OR.XW,EO.XE>  GO  TO  13 
DO  12  J=l<  TMAX 

12  Y  (  1 1  J  )  =Y  I  (  IWtJ)-MYI  (  IEtJ)-YI  <  IWf  J)  )  *  (  X  (  I  ) -XW  )  /  (  XE-XW  ) 

GO  TO  15 

13  IW=IW+1 
IE»IW+1 
XW=XB( IW) 

XE=XB( IE) 

GO  TO  11 

15  CONTINUE 

DO  20  1=1 i IMAX 
20  WRITE<  8 1 200 )  X<I) 

DO  30  J=  1 1  JMAX 
30  WRITE (8i 200 )  Z(J) 

DO  40  I»lf IMAX 
DO  40  J=ltJMAX 
40  WRITE ( 8t  230 )  Y(ItJ) 

100  FORMAT  <  3F 1 0 . 5 i  11X»2F13,5) 

110  FORMAT  (  1 0X  i  2F 10.  5 1  HXi2F10.5) 

200  FORMAT  < 1  PEI 3. 6 ) 

END 

SUBROUTINE  NAMIN 

C  ****•»#*#•»*#•**#*****##**********■»**♦***•»#*#*#*#**#*#****#*### 


06100  C  READ  INPUT  DATA  FILE 
06200  C  *******•*•*■*■*•*•*■**•**■-*■*■* 

06300  C  DESCRIPTION  OF  ’ PRMESHIN’  INPUT  DATA: 

06400  C  IMAX. . .NUMBER  OF  POINTS  IN  STREAMWISE  DIRECTION:  OUTPUT 

06500  C  JMAX. . .NUMBER  OF  POINTS  IN  VERTICAL  DIRECTION:  OUTPUT 

06600  C  NXTL. .. NUMBER  OF  POINTS  IN  STREAMWISE  DIRECTION:  INPUT 

06700  C  NZTI _ NUMBER  OF  POINTS  IN  VERTICAL  DIRECTION:  INPUT 

06800  C  X( I )... SPECIFIED  STREAMWISE  COORDINATE  ARRAY:  OUTPUT 

06900  C  Z (J) ...  SPECIFIED  VERTICAL  COORDINATE  ARRAY:  OUTPUT 

07000  C  ******************************-*****-»******************#***-** 
07100  COMMON/DATA/  X(41)»Z(20) 

07200  COMMON/NDATA/  N(4) 

07300  CALL  OPEN  (  7,  ’  PRMESHIN  M) 

07400  CALL  OPEN  <8*  ’  MESHPR  M) 

07500  CALL  OPEN! 9,’ MESH  ’,1) 

07600  DO  1  1=1*4 

07700  N ( I )  =  0 

07800  1  CALL  NAMLST  <7*2*N(I)*  DUMMY ) 

07900  IMAX=N(1) 

08000  JMAX=N!2) 

08100  DO  2  1  =  1  * IMAX 

08200  2  CALL  NAMLST (7*2* 1  *  X  !  I  )  ) 

08300  DO  3  J=1*JMAX 

08400  3  CALL  NAMLST (7*2* 1*Z(J) ) 

08500  RETURN 

08600  END 

08700  SUBROUTINE  NAMLST! 1 1  *  10, N* X ) 

08800  C  **#***#*****#****»***********♦****#*♦**##****#***»*****«#*** 
08900  C  READ  vaRIaBlE  NAME  AND  VALUE 

09000  C  #*****************■*■*■****************•*•*■*•**■*■■*■*•**•**■■**•*•*****•**** 
09100  DIMENSION  A<3> 

09200  IF ( N. NE. 0 )  GO  TO  1 

09300  READ! 11,2)  A!1)*A!2)»A!3)*N 

09400  WRITE! 10, 2)  A!1),A!2)*A!3),N 

09500  RETURN 

09600  1  READ ! 1 1  *  3 )  A !  1  )  ,  A  !  2  >  *  A  !  3 ) ,  X 

09700  WRITE ( 10,3)  A ! 1 ) , A ! 2 ) *  A ! 3 ) *  X 

09800  RETURN 

09900  2  FORMAT! IX, 3A4* 14) 

10000  3  FORMAT! IX, 3A4, 1PE13. 6) 

10100  END 


300 

310 

320 

330 

340 

350 

360 

370 

380 

390 

400 

410 

420 

430 

440 

450 

460 

470 

480 

490 

500 

510 

520 

530 

540 

550 

560 

570 

580 

590 


>  REM  ******************.**.*.**.*.*., 

>  REM  **********  PLOTTING  ROUT] 

>  REM  **********  DISK  FILE  *1" 

REM  **********  PROGF 

REM  **********************<Mhh, 
DIM  X(41 ) » Y<20> » Z<20> 

NZ=70 

NY=20 

NL=9 

IC=0 

OPEN  "I",  1,  ■ MESHPR" 

INPUT  "IMAX  =  "  !  IM 
INPUT  “ JMAX  =  •  ?  JM 
FOR  1=1  TO  IM 
INPUT#li  X<I> 

NEXT  I 

FOR  J=1  TO  JM 
INPUT#1,  Z(J) 

NEXT  J 
IC=IC+1 

IF  I C > I M  THEN  1600 
FOR  J=1  TO  JM 
INPUT#1 »  Y(J> 

Y<  J) =— Y ( J ) 

NEXT  J 
GOSUB  620 
CLS 

PRINT  “I  =" ;  IC;  *  X<I)  =■ 

on  r  mt ■  noTMT 

•  «  \  A  I  *  •  -  I  I  \  J.  <  1  I 

PRINT  " Z  M I N  =  “5  ZL?  -  ZMAX 
PRINT 

PRINT  " YMIN  =  "?  YL?  •  YMAX  < 
PRINTJPRINT 

PRINT  "SPECIFY  PLOTTING  RANGE* 
PRINT 

IF  IC=1  THEN  430 
PRINT  "SAME  RANGE  ?" 

Q*= INKEY* 

IF  0*=""  THEN  380 

IF  LEFT* ( 0*i 1 ) <>«Y"  THEN  430 

GOSUB  1550 

GOTO  510 

INPUT  "ZMIN  =« ;  ZL 
INPUT  "ZMAX  =  "?  ZU 
INPUT  "YMIN  =  •?  YL 
INPUT  "YMAX  =  "?  YU 
GOSUB  1490 
PRINT 

INPUT  "ANY  ERRORS" !  0* 

IF  LEFT* ( 0*i 1 )="Y"  THEN  340 
CLS 

GOSUB  770 
GOSUB  1060 
GOSUB  1130 
A*= INKEY* 

IF  A*= " "  THEN  550 
IF  A*="C"  THEN  260 
IF  A*="N"  THEN  200 
END 


PLOTTING  ROUTINE  FOR  SHIP  SECTIONS  ********** 
DISK  FILE  ’MESHPR’  CREATED  BY  ********** 
PROGRAM  PRMESH  ********** 


X<I)  =" 5  X(IC> 
ZMAX  =  “ 5  ZU 


YMAX  = 


:  0* 

THEN  340 
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600  REM  **********  FIND  MIN  &  MAX  VALUES  ********** 

610  REM  -  MIN  VALUE  - 

620  ZL-1.E20 

630  YL-1.E20 

640  FOR  1  =  1  TO  JM 

650  IF  Z(IXZL  THEN  ZL=Z<I) 

660  IF  Y(IXYL  THEN  YL=Y<I) 

670  NEXT  I 

680  REM  -  MAX  VALUE  - 

690  ZU=ZL 

700  YU=YL 

710  FOR  1=1  TO  JM 

720  IF  Z<I)>ZU  THEN  ZU=Z(I) 

730  IF  Y < I ) > YU  THEN  YU=Y(I) 

740  NEXT  I 
750  RETURN 

760  REM  **********  DRAW  AXES  AND  TIC  MARKS  ********** 
770  ZZ=ZU— ZL 
780  YY=YU— YL 

790  REM  -  Y-AXIS  TIC  MARKS  - 

800  IY-/.=  1 

810  FOR  1=1  TO  5 

820  PR  I  NT  a)  (  I YX  *  NL— 2  )  i  CHR*<95); 

830  PRINTS*  (  IY*/.i  NL— 1  )  i  CHR*<95> 

840  IY7.=  IY7.+5 
850  NEXT  I 

860  REM  -  DRAW  THE  Y-AXIS  - 

870  IY7.=2 

1  1  =NY+i 

890  FOR  I = I Y%  TO  II 

900  PRINTS ( I »  NL ) f  CHR*<156) 

910  NEXT  I 

920  REM  -  Z-AXIS  TIC  MARKS  - 

930  1 1 =NZ+NL 
940  IS=NZ/5 

950  FOR  I =NL  TO  II  STEP  IS 
960  PRINTS  <  NY+2i  I  )  »  CHR*<156)? 

970  NEXT  I 

980  REM  -  DRAW  THE  Z-AXIS  - 

990  I I=NZ+NL— 1 
1000  NN=NL+1 
1010  FOR  I =NN  TO  II 
1020  PRINTS ( NY+ 1  *  I ) i  CHR*<95) 

1030  NEXT  I 
1040  RETURN 

1050  REM  **********  LABEL  AXES  ********** 

1060  PRINT S ( 0*  0 )  *  YU? 

1070  PRINTS ( NY+ 1  *  0 )  i  YL ? 

1080  PRINTS ( NY+3»  NL— 2 )  »  ZL ? 

1090  PRINTS  <  NY+3*  NZ+NL— 6 )  i  ZU? 

1100  PRINTS ( 0i NZ /2+NL— 4 ) *  "X  =  ■  ?  X(IC)I 
1110  RETURN 

1120  REM  **********  PLOT  THE  POINTS  ********** 

1130  ID=0 

1140  FOR  1=1  TO  JM 

1150  C»NL+INT<NZ*< Z ( I ) -ZL ) /Z Z+. 5 ) 

1160  R=NY+ 1- INT  <  NY*  <  Y  < I ) -YL ) /YY ) 

1170  IF  C> (NZ+NL)  THEN  1270 
1180  IF  C<NL  THEN  1270 
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1190  IF  RXNY+t)  THEN  1270 
1200  IF  R<  1  THEN  1270 
1210  GOSUB  1360 
1220  IF  GC*=160  THEN  1240 
1230  PRINTS <  Ri C ) i  CHR*<GC>? 

1240  PRINTS  <  R+l »  C ) »  CHR$  <  26 ) +CHR*  <  GC ) ; 

1250  PRINT  CHR*<25>; 

1260  GOTO  1280 

1270  ID=ID+1 

1280  NEXT  I 

1290  PRINTS ( 7i 1 ) t  ID? 

1300  PRINTS ( 8i 0 ) t  "POINTS"? 

1310  PRINTS ( 9*  1 ) f  "OFF"; 

1320  PRINTS < 10i0) i  "SCALE"; 

1330  PRINT  CHR*<02); 

1340  RETURN 

1350  REM  **********  INTERPOLATE  WITHIN  BOX  ********** 

1360  RR=R— ( NY+ 1— NY*  <  Y ( I ) — YL ) /YY ) 

1370  GC—154 

1380  IF  RR< . 2  THEN  RETURN 
1390  GC= 155 

1400  IF  RR< . 4  THEN  RETURN 
1410  R=R— 1 
1420  GC— 160 

1430  IF  RR< . 6  THEN  RETURN 
1440  GC=152 

1450  IF  R  R  <  •  8  THEN  RETURN 
1460  GC=1 53 
1470  RETURN 

1480  REM  **********  SAVE  PLOTTING  RANGES  ********** 

1490  XL—ZL 
1500  XU^ZU 
1510  WL=YL 
1520  WU=YU 
1530  RETURN 

1540  REM  **********  RETRIEVE  PLOTTING  RANGES  ********** 

1550  ZL=XL 

1560  ZU=XU 

1570  YL=WL 

1580  YU=WU 

1590  RETURN 

1600  END 
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PROGRAM  SHIP 

C  ********************************************************************* 
C  THIS  PROGRAM  USES  LINEAR  SHIP  THEORY  TO  COMPUTE  THE  WAVE  DRAG 
C  AND  SURFACE  VELOCITY  COMPONENTS  FOR  AN  ARBITRARY  VESSEL. 

C  ********************************************************************* 
COMMON/COORD/  XB(41 ) « YB(41 «  10) «  ZB( 10) « FB(40*9> 

COMMON/DATA/  FRi  IMAX* JMAX* MINT2* MINT3* NZOT* IMAXMi  JMAXM 
COMMON/QUADRA/ AM <  80  >  i RSAV (80)* TSAV  <  80  > i WM 
COMMON/ VELO/  U ( 40 t 9 ) f W < 40, 9 ) 

CALL  NAM IN 
FR2R=1 . /FR**2 
DO  1  1  =  1 1 IMAX 
READ ( 7 1 100 )  XB  <  I  ) 

1  XB ( I )=FR2R*XB< I ) 

DO  2  J=1 ?  JMAX 
READ ( 7 1  100 )  ZB  <  J ) 

2  ZB< J)=FR2R*ZB< J) 

DO  3  1  =  1, IMAX 

DO  3  J=1,JMAX 
READ ( 7 i 100)  YB ( I* J) 

3  YB  < I ,  J ) =FR2R*YB ( I?  J) 

CALL  GEOM< IMAX i JMAX ) 

RWAVE=0. 

DO  4  1=1 t IMAXM 
DO  4  J=1,JMAXM 
U<I,J)=0. 

4  W(I,J>=0. 

CALL  VELOC1 
Wn I TE ( i »  3 10 ) 

CALL  VEL0C2 
WRITE  <  1  >  520 ) 

CALL  VEL0C3 
WRITE  < 1  *  530 ) 

DO  10  1  =  1  *  IMAXM 
CALL  DRAG ( DRWAVE, I ) 

RWAVE=RWAVE+D RWAVE 

IF  <  NZOT. EQ.  0 )  CALL  OUTPUT ( I ) 

WRITE (2* 400)  DRWAVE, RWAVE 
10  CONTINUE 

WRITE (2i 200)  RWAVE 
WRITEClf 200)  RWAVE 
100  FORMAT ( E13- 6 ) 

200  FORMAT (////// i 1 X 1 8HC0MPUTED  WAVE  DRAG/ / /4X7HRWAVE  =,1PE12.4> 

400  FORMAT <//, IX t • INCREMENTAL  WAVE  DRAG  ON  THIS  SECTION  IS  ?,1PE13.4/ 

*  ,  1  X  i  ’ TOTAL  DRAG  COMPUTED  UP  TO  THIS  SECTION  IS’,E13.4) 

510  FORMAT ( 1X» ’ VELOC 1  COMPUTATIONS  COMPLETE* ) 

520  FORMAT ( 1 X*  * VEL0C2  COMPUTATIONS  COMPLETE’) 

530  FORMAT (IX*’ VEL0C3  COMPUTATIONS  COMPLETE’) 

END 

FUNCTION  ACOS  <  X ) 

C  ************** ******************************************************* 
C  THIS  SUBROUTINE  COMPUTES  THE  ARC  COSINE 

C  ********************************************************************* 
DATA  PH/ 1 . 3707963268/ 

ACOS*PH— AS I N <  X ) 

RETURN 

END 

FUNCTION  AS IN ( X  > 

C  ********************************************************************* 
C  THIS  SUBROUTINE  COMPUTES  THE  ARC  SINE 
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C  ***********************■**#****•***************•*****#***********#*****’ 
DATA  A0i A 1 i A2/ 1 . 57079630501 2145988016, . 0889789874/ 

DATA  A3i A4i A5/-. 0501743046* . 030891 8810* 017088 1256/ 

DATA  A6, A7»  PH/ . 0066700901 001 262491 1,1. 3707963268/ 

Y*ABS<X) 

AS  I N-PH-SQRT  (  1 .  — Y  )*<  A0+Y*  (  A 1 +Y*  ( A2+Y*  (  A3+ Y*  ( A4+Y*  <  A5 
*  +Y*(A6+Y*A7> ) ) ) > ) ) 

IF <  X  • LT • 0. )  AS IN*— AS IN 

RETURN 

END 

SUBROUTINE  DRAG ( DRWAVE, I ) 

C  **************«^********#***************#***#*************#********** 
C  THIS  SUBROUTINE  COMPUTES  THE  WAVE  DRAG 

C  ******************************************************************** 
COMMON/COORD/  XB(41 )i YB(41i 1 0 ) i ZB ( 10 ) t FB ( 40i 9 ) 

COMMON/DATA/  FRi IMAX i JMAX i MINT2* MINT3i NZOTi IMAXMi JMAXM 
COMMON/VELO/  U ( 40 t 9 ) i W ( 40* 9 ) 

DRWAVE=0. 

DX«XB< 1+1 >-XB< I ) 

DO  10  J=2, JMAX 
DZ«ZB< J-l )-ZB< J) 

10  DRWAVE=DRWAVE+U ( 1 1 J- 1 ) *FB ( I  *  J- 1 ) *DX*DZ 
DRWAVE— 4 . *DRWAVE*FR**4 
RETURN 
END 

SUBROUTINE  GEOM < IMAX , JMAX > 

C  ************************************************************************** 
C  THIS  SUBROUTINE  COMPUTES  CELL  SLOPE 

COMMON/ COORD/  XB< 41 ) , YBC41 , 10) , ZB< 10) iFB(40i 9) 

DO  10  I =2* IMAX 
DX«XB< I ) — XB ( I— 1 ) 

DO  10  J=2t  JMAX 

10  FB ( I  —  1 1 J- 1 )  * .  5  *  ( YB  < 1 1 J)+YB< 1 1 J-l )-YB< 1-1 i J)-YB( 1-1 1 J-l ) > /DX 

RETURN 
END 

SUBROUTINE  INTK0 ( G*  Z i RMU ) 

C  *****************************************************-*•#■*■*•**■*■*•***•*■*** 
C  THIS  SUBROUTINE  APPROXIMATES  THE  BESSEL  FUNCTION  INTEGRAL  WITH  A 
C  GROUP  OF  CLOSED  FORM  ANALYTICAL  FUNCTIONS. 

C  ******************************************************************** 
ZMU=Z*RMU 

IF <  ZMU. LE. 0 )  GO  TO  20 

IF ( ZMU. LE. .5)  G=ZMU* ( 1 . — ALOG <  ZMU ) ) 

IF <  ZMU. GT.  . 5 )  G*1.2*S0RT< ZMU) 

ZMU2*ZMU*RMU 

IF  <  Z . LT. 3.  )  GO  TO  40 

X 1*3. /Z 

GMAX*1 .  +.  571* < l.-XI >*(.37+.63*EXP(-10.*XI )  ) 

IF(GMAX.LT.G)  GO  TO  10 
RETURN 


10  G*GMAX 

IF ( Z . GE. 10. )  GO  TO  80 
GO  TO  60 

C - MU* Z  IS  ZERO - 

20  G*0. 

RETURN 

C - Z  BETWEEN  0  AND  3  - 

40  RMUT  * 1 . ♦ 1 • 114*Z**<-. 6666667) 

I F  <  RMU • GE • RMUT )  GO  TO  35 


12100  I F  (  Z  . GE . 1 >  GO  TO  45 

12200  ZP25=Z**.25 

12300  GM=.83*ZP25 

1 2400  RN= . 62+ 1 . 60*  Z  P25 

12500  GO  TO  50 

12600  45  GM=.83+<Z-1. )*(. 12-.02*(Z-1. ) ) 

12700  RN=2.+.3*Z 

12000  50  G2=GM/(1.+.25*ZMU2**RN> 

12900  GO  TO  100 

13000  55  G=SQRT< 1. 5707963/ ZMU>*EX P(-ZMU) /( RMU- 1. ) 

13100  RETURN 

13200  C  -  Z  BETWEEN  3  AND  10  - 

13300  60  RMUT=l.+l. 114* Z**<-. 6666667) 

13400  IF<  RMU. GE. RMUT )  GO  TO  65 

13500  XI=.l*(Z-3.) 

13600  A=1 . 30+X I* ( . 76—. 54* X I ) 

13700  RM=. 68+X I*  < . 53—. 41 *X I ) 

13800  G2=A*EXP ( — RM*ZMU2 ) 

13900  GO  TO  100 

14000  65  G=SQRT ( 1 . 5707963/ZMU ) *EXP (  —  ZMU ) / <  RMU— 1 .  ) 

14100  RETURN 

14200  C  -  Z  GREATER  THAN  10  - 

14300  80  G2— 0. 

14400  IF <  RMU. LT. 1 .  )  G2= < 1 . 571+ASIN < RMU ) ) *EXP < -ZMU2 > 

14500  *  /SORT ( 1 . — RMU*RMU ) 

14600  100  IF ( G2. LT . G  >  G=G2 

14700  RETURN 

14800  END 

14900  SUBROUIlNE  NAM IN 


15000  C  ************************************************************ 
15100  C  READ  INPUT  DATA  FILE 
15200  C  ******************** 

15300  C  DESCRIPTION  OF  ’SHIPIN'  INPUT  DATAi 


15400  C  FR . FROUDE  NUMBER  *  U/SQRT<G*L> 

15500  C  IMAX. . .NUMBER  OF  MESH  POINTS  IN  STREAMWISE  DIRECTION 

15600  C  JMAX. .. NUMBER  OF  MESH  POINTS  IN  VERTICAL  DIRECTION 

15700  C  MINT2. .NUMBER  OF  POINTS  IN  QUADRATURE. .. INFINITE  INTEGRAL 

15800  C  MINT3. .NUMBER  OF  POINTS  IN  QUADRATURE . FINITE  INTEGRAL 

15900  C  NZOT...0  TO  PRINT  CELL  GEOMETRY  AND  VELOCITY  COMPONENTS 

16000  C  ...OTHERWISE  PRINT  ONLY  THE  WAVE  DRAG 

16100  C  ************************************************************ 
16200  COMMON/DATA/  FR, N ( 5 > , IMAXM, JMAXM 

16300  CALL  OPEN (7, ' MESHPR  ’,1) 

16400  CALL  OPENO,  ’SHIPIN  ’,1) 

16500  CALL  OPEN (9, ’  UEWE  ’,1) 

16600  CALL  NAMLST (8, 2, 1 , FR ) 

16700  DO  1  1=1,5 

16800  N( I )  =0 

16900  1  CALL  NAMLST (8,2,N(I), DUMMY ) 

17000  IMAXM=N ( 1 ) — 1 

17100  JMAXM=N ( 2 ) — 1 

17200  RETURN 

17300  END 

17400  SUBROUTINE  NAMLST ( 1 1 , 10, N, X ) 

1 7500  C  ************************************************************ 
17600  C  READ  VARIABLE  NAME  AND  VALUE 

1 7700  C  ************************************************************ 
17800  DIMENSION  A(3> 

17900  IF ( N. NE. 0 )  GO  TO  1 

18000  READ (11,2)  A( 1 ) , A(2) , A(3) ,N 

-Hi- 


18100 

18200 

18300 

18400 

18300 

18600 

18700 

18800 

18900 

19000 

19100 

19200 

19300 

19400 

19500 

19600 

19700 

19800 

19900 

20000 

20100 

20200 

20300 

20400 

20300 

20600 

20700 

20800 

209MV9 

21000 

21100 

21200 

21300 

21400 

21300 

21600 

21700 

21800 

21900 

22000 

22100 

22200 

22300 

22400 

22300 

22600 

22700 

22800 

22900 

23000 

23100 

23200 

23300 

23400 

23300 

23600 

23700 

23800 

23900 

24000 


WRITE  <  I0»  2 )  A(l),A(2)iA(3)iN  < 

RETURN 

1  READ  <IIi3>  A ( 1 ) « A ( 2 ) v A ( 3 > v X 
WRITE  < IOi 3 )  A( 1 )iA(2)iA(3)iX 
RETURN 

2  FORMAT (  1  X »  3A4i  14) 

3  FORMAT  <  1  X  i  3A4»  1PEK3.6) 

END 

SUBROUTINE  OUTPUT<I) 

C  ********************************************************************* 
C  THIS  SUBROUTINE  PRINTS  SURFACE  VELOCITY  COMPONENTS  AND  CREATES  DISK 
C  FILE  ’ UEWE*  FOR  PROGRAM  VELOC 

C  ********************************************************************* 
COMMON/COORD/  XBC41 > f YBC41 *  10) «  ZB( 10) «FB(40v9> 

COMMON/DATA/  FR, IMAX  »  JMAX t  MINT2i  MINT3* NZOT*  IMAXMt  JMAXM 
COMMON/VELO/  U(40i 9) i W<40,9) 

FR2=FR*FR 

XI».5*FR2*(XB< I )+XB< 1  +  1 ) ) 

WRITE  <  2i 100 )  XI 
DO  1  J=1 »  JMAXM 
JJ-JMAX-J 

ZETA«.5*FR2*(ZB< JJ)+ZB< JJ+1 ) > 

ETA=.5*FR2*<YB< I iJJ )+YB( I* JJ+1 ) ) 

UP=1.+U< If JJ) 

VP=FB ( I f  J J ) 

WP=W( I  *  JJ) 

WRITE  <  2t  200 )  Jf ZETAf UPi VPf WP 
Q-SQRT ( UP*UP+VP*VP+WP*WP ) 

UP=UF/Q 

VP-VP/Q 

WP=WP/Q 

1  WR I TE ( 9i  1000)  XI i ETA f  ZETA*  UP*  VP» WP?Q 

100  FORMAT  (  / 1  1  X38HC0MPUTED  VELOCITY  COMPONENTS  ALONG  X  *ilPE12.4// 

*  4X1 HJf 7X4HZ< J) ,  1 1 X4HU ( J  >  » 11X4HV(J)* 11X4HW(J)/) 

200  FORMAT ( I5» 1P4E15. 4) 

1000  FORMAT ( 3F 1 0 .  5  i 4F 10.6) 

RETURN 

END 

SUBROUTINE  OUAD(M) 

C  *********************************************************.**.**.*.**.*.**.* 
C  THIS  SUBROUTINE  COMPUTES  QUADRATURE  WEIGHTS  AND  ABSCISSAS 
C  #*#***#***#****#*******#************»#*#******#*******•**’'**«’***-*•»*•»** 
COMMON/DATA/  FRi  IMAX i JMAX  »  MINT2*  MINT3i NZOT *  IMAXM* JMAXM 
CO*  MON /QUADRA /AM  (  80  )  i  RSAV  (  00  )  i  TSAV  (  80  )  t  WM 
DATA  PI/3.  14159265/*  CON/ . 2026423/ 

C  EQUAL  WEIGHT  CHEBYSHEV  INTEGRATION  FOR  DU3  AND  DW3 
IF  <  M. EQ. 2 )  GO  TO  15 
COEF*P I / ( 4 . *FLOAT ( M I NT3 ) ) 

WM*2. *CON*COEF 
DO  10  1*1 i MINT3 
RI=FLOAT( I ) 

AM< I )*COS( ( 2. *RI— 1 • )*COEF) 

RSAV  < I ) =SQRT  < 1 . -AM ( I ) **2 ) 

10  TSAV  < I ) *ACOS ( AM ( I )  )+PI 
RETURN 

C  EQUAL  WEIGHT  CHEBYSHEV  INTEGRATION  FOR  DU2  AND  DW2 
15  COEF*PI / ( 4 ■ *FLOAT  <  MINT2 ) ) 

WM*2.*C0N*C0EF 
DO  20  I»1 ?  MINT2 
R I -FLOAT < I ) 
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AM< I ) “COS ( (2.*RI-1. >*COEF> 

EM2«AM( I )**2 

RSAV ( I ) -SORT  < 1 . -EM2 ) /AM ( I ) 

ARGLOG**!  .  /AM  (  I  )  -RSAV  (  I  ) 

20  TSAV ( I ) “ALOG <  ARGLOG ) 

RETURN 

END 

SUBROUTINE  VELOC1 

C  ****************************************************************************** 
C  THIS  SUBROUTINE  COMPUTES  THE  CONTRIBUTION  OF  THE  SOURCE-LIKE  TERMS 
C  TO  THE  SURFACE  VELOCITY  COMPONENTS 

C  *********************************************■**■******#******■**-*.*.*#■**# 
COMMON/COORD/  XB ( 41 >  *  YB ( 41  *  10 >  *  ZB ( 10 >  * FB ( 40* 9 > 

COMMON/DATA/  FR, IMAX * JMAX*  MINT2*  MINT3?  NZOT 9  IMAXM*  JMAXM 
COMMON/VELO/  U(40*  9) 1 W(40* 9) 

DATA  PI/3. 14159265/ 

C  -  INITIALIZATION  - 

DO  40  1=1 * IMAXM 
XI*=.  5*(  XB<  I  >+XB(  I  +  l)  ) 

DO  40  J=l* JMAXM 
ZETA**.  5*  (  ZB(J)+ZB(J+1  )  ) 

UINT=0. 

WINT=0. 

DO  30  11  =  1* IMAXM 
DO  30  JJ=1* JMAXM 
XP=XB (II) —X  I 
XM=*XB  <  I  1  +  1  )  —  X  I 
ZP=- <  ZETA+ZB ( J  J+ 1 ) ) 

ZM«-<ZETA+ZB< JJ) ) 

Z 1 P=ZB ( JJ+  1 ) — ZETA 
Z 1M=ZB  (JJ)— ZETA 
RA=SQRT<  XP*XP+ZP*ZP) 

RB*=SQRT  <  XM*XM+ZP*ZP ) 

RC»SQRT<  XM*XM+ZM*ZM) 

RD=SQRT( XP*XP+ZM*ZM) 

R1 ASSORT ( XP*XP+Z1P*Z1P> 

R1B=SQRT ( XM*XM+Z 1 P*Z 1 P ) 

R1 C“SQRT ( XM*XM+Z 1M*Z 1M ) 

RlD=SORT ( XP*XP+Z 1M*Z 1M ) 

UINT*UINT+ALOG<  < R1 A-Z 1 P ) * <  R1 C-Z 1 M ) * <  RD+ZM ) * ( RB+ZP > 

*  /( ( RIB— Z 1 P ) * ( RID— Z 1M ) * <  RC+ZM )  #  ( RA+ZP )  >  >*FB< II?  JJ ) 

W I NT «W I NT +ALOG  <  <  R 1 A-XP ) * < R 1 C-XM ) * < RD-XP ) * ( RB-XM ) 

*  /( ( R1 D— XP ) *  <  RIB—  XM ) *  <  RC— XM ) *  <  RA— XP )  )  >*FB< II* JJ) 

30  CONTINUE 

U(I* J)«U( It J)+.3*UINT/PI 
W< It J)*W( It J)+.3*WINT/PI 
40  CONTINUE 
RETURN 
END 

SUBROUTINE  VEL0C2 

C  **#■***■»*##**•»##**«■***♦***#*#*♦****-»-»***♦**********#********#****•»***♦ 
C  THIS  SUBROUTINE  COMPUTES  THE  CONTRIBUTION  OF  THE  INFINITE  INTEGRAL 
C  TO  THE  SURFACE  VELOCITY  COMPONENTS 

C  *■»#*******#*-#***#*#**♦**  ************  *•»**#***  ***MMMMMf*  *****•*■*##*■*•****« 

COMMON/COORD/  XB < 41 ) t YB < 41 1 10 ) t  ZB ( 10 > t FB ( 40 1 9 > 

COMMON/DATA/  FR* IMAX * JMAX * MXNT2* MINT3* NZOT* IMAXM* JMAXM 
COMMON/ I NTG/  TC0S<41 ) *  TSIN(41 ) *  TEXP(9*  10 >  *  TBES(9*  10 > 
COMMON/OUADRA/AM ( 80 ) * RSAV (00) *  TSAV (80) *  WM 
COMMON/VELO/  U ( 40*  9 ) *  W ( 40*  9 ) 

DATA  PI/3.14139263/ 
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C  - - INITIALIZATION - 

CALL  QUAD ( 2 ) 

DO  100  M=1 , MINT2 
EM*AM<M) 

EM2=EM*EM 
DO  10  J=l, JMAXM 
ZETA*. 5*  <  ZB(J)+ZB<J+1 ) ) 

DO  10  JJ=1,JMAX 
ZED—— ( ZB  <  J J ) +ZETA ) 

EMR=1 . /EM 

CALL  INTK0 (6i ZED,  EMR ) 

TBES ( J, J J ) =G*RSAV ( M ) 

10  TEXP  <  J,  JJ ) =EXP  < -ZED/EM2 ) 

C  -  COMPUTE  FOURIER  INTEGRAL  - 

DO  40  1=1, IMAXM 
XI=.5*(XBCI  )  +  XB(  14-1)  ) 

DO  20  1 1  =  1 ,  IMAX 
T COS ( 1 1 )=COS( <  XB< 1 1 >-XI >/EM) 

20  TSIN< 1 1 ) =SIN ( <XB( II )  — XI >/EM) 

DO  40  J=l, JMAXM 
UINT=0. 

WINT =0. 

DO  30  11=1, IMAXM 
COSXOM=T COS  < II )-TCOS< II  +  l ) 

SINXOM=TSIN( I  I ) -TSIN  < 1 1  +  1 ) 

DO  30  J J= 1 i JMAXM 

C=<TEXP< J, JJ+1 ) — TEXP ( J ,  J  J  )  >*FB< I  I , J J ) 

D=  <  TBES  ( J,  JJ+ 1  >  —TBES  (  J ,  JJ )  ) *FB C 1 1 , J J > 

U I  NT =U  I  NT C*  (  PI*5INX0m-7SAv  i  h  ;  *C06*0M ;  —  U*cubXon 
30  WINT=WINT+  <  C* ( PI*COSXOM-TSAV  <  M ) *S INXOM  > -D*S INXOM > /EM 

U<  I,  J)=U(  I,  J)+WM*UINT 
W( I, J)=W( I, J)+WM*WINT 
40  CONTINUE 
100  CONTINUE 
RETURN 
END 

SUBROUTINE  VEL0C3 

C  ******************************************************************** 
C  THIS  SUBROUTINE  COMPUTES  THE  CONTRIBUTION  OF  THE  FINITE  INTEGRAL 
C  TO  THE  SURFACE  VELOCITY  COMPONENTS 

C  ******************************************************************** 
COMMON/ COORD/  XB<41 ) , YB<41 , 10) , ZB( 10) , FB<40, 9) 

COMMON/DATA/  FR,  IMAX , JMAX , MINT2, MINT3, NZOT  *  IMAXM, JMAXM 
COMMON/ I NTG/  TC0S<41 ) , TSIN <41 ) , TEXP (9, 10) , TBES (9, 10) 
COMMON/QUADRA/AM ( 80 ) , RSAV ( 80 ) , TSAV ( 80 ) , WM 
COMMON/ VELO/  U < 40 , 9 ) , W ( 40 *  9 ) 

C  -  INITIALIZATION  - 

CALL  QUAD ( 3 ) 

DO  100  M= 1  ,  MINT3 
EM=AM ( M ) 

EM2=EM*EM 

DO  10  J=l, JMAXM 

ZETA=, 5*( ZB<J)+ZS(J+1 ) ) 

DO  10  JJ= 1 , JMAX 
ZED*— ( ZB ( JJ ) +ZETA ) 

CALL  INTK0 ( G, ZED, EM ) 

1 0  TEXP  <  J , J J ) =TSAV ( M ) *EXP  < -ZED*EM2 ) +RSAV ( M ) *G 

C  -  COMPUTE  FOURIER  INTEGRAL  - 

DO  40  1=1, IMAXM 

XI  =  .5*<XB(  I  )+XB<  1 4*  1  )  ) 


36100 

DO  20  I 1*1 v IMAX 

36200 

T  COS (II) —COS  < <XB< II >-XI >*EM> 

36300 

20 

TSIN( 1 1 )=SIN< <  XB < II >-XI >*EM> 

36400 

DO  40  J=1»JMAXM 

36300 

UINT=0. 

36600 

WINT=0. 

36700 

DO  30  11=1 » IMAXM 

36800 

COSXMU=*T COS  <  1 1  )  — T COS  (  1 1  + 1  ) 

36900 

SINXMU=TSIN< II >-TSIN< II+l ) 

37000 

DO  30  JJ=1 1  JMAXM 

37100 

C=<TEXP< J* JJ+l ) — TEXP ( Ji JJ ) >*FB( II* JJ) 

37200 

UINT=UINT-C*COSXMU/EM 

37300 

30 

WINT=WINT+C*SINXMU 

37400 

U< I* J)=U< I* J)+WM*UINT 

37300 

W< I* J)=W( I* J)+WM*WINT 

37600 

40 

CONTINUE 

37700 

100 

CONTINUE 

37800 

RETURN 

37900 

END 
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PROGRAM  SHPMSH 

C  *****  »##**•»*#***#*#**#***#*#*****  »***#**##*######**  » #**«**«•»«-***#*** 
C  THIS  PROGRAM  CALCULATES  METRIC  AND  CURVATURE  PARAMETERS  FOR  ARBITRAR 
C  SHIP  HULLS  AND  WRITES  A  DISK  FILE  FOR  USE  BY  PROGRAM  EDDY3 
C  ***************************************************  ****************** 
COMMON/COORD/  Z < 40 ) 

COMMON/ COORDB/XB ( 40  >  t  YB ( 40 ) i ZB ( 40 ) 

COMMON/ INTRP/  X  I ( 40 > » Z I ( 40 > 

COMMON/GEOMP/  H 1 ( 40 ) ?  H2 ( 40 ) i  THETA ( 40i 40 ) *  CK 1 ( 40 ) t CK2 <  40 ) » 

*  CK12 ( 40 ) ? CK21 (40) , PSI (40) 

DYBDZ ( 40 ) i D2YBDZ <  40 ) * DZBDZ ( 40  >  t  D2ZBDZ  <  40 ) 

DYBDX ( 40 ) i D2YBDX ( 40  >  i  DZBDX  <  40 ) *  D2ZBDX ( 40  > 

DTHDX <  40 ) »  DTHDZ < 40 ) * THET ( 40 ) *  LT ( 30 ) 

PSIZi XSTARTt  XENDi DX < 10 >  v NX < 10 > 

ZSTARTt  ZEND »  DZ ( 10) f NZ < 10) 

IUlt  IU2i  IU3i  IU4»  IU5i I MAX  $  JMAX  t  NXTLi  NZTLi  JFZ iLZ ( 30 ) 
JFX  i  LX ( 30 ) 

ISYMt I SET 

i  1  ) 


COMMON/TEMPZ / 
COMMON/ TEMP X/ 
COMMON/DTHET/ 
COMMON/ XDATA/ 
COMMON/ ZDATA/ 
COMMON/NDATA/ 
COMMON/MDATA/ 
COMMON/LDATA/ 
CALL 
CALL 
CALL 
CALL 
CALL 


OPEN <6, ’ NAMLIST 
NAM  IN 

OPEN  < IUli  ’MESH  ’ i  1 ) 
OPEN ( I U2 i ’SCRATCH1 
OPEN ( IU3i  ’ SC RAT  CH2 


CALL  OPEN ( IU4i ’ GEOMETRY 
SET  UP  COMPUTATIONAL  GRID 
CALL  GRID 
JMAXM=JMAX-1 
DO  1  I = 1 i I MAX 

X  I 


i  12) 

i2S) 
i  32) 


03000  C 
03100 
03200 
03300 
03400  2 


READ  SECTION  COORDINATES 

READ (IUli  1000)  XB ( I ) i YB (  1 ) »  ZB( 1 ) » YB < 2 ) i ZB ( 

DO  2  J*=3i  JMAXMi  2 

JP=J+1 

READ ( IU1 » 1100)  YB( u) i ZB( J) i YB( JP) i ZB( JP) 


) 
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C  COMPUTE  AND  RESCALE  LOCAL  ARCLENGTH 
CALL  ARC ( DYDX 1 ) 

C  SMOOTH  SECTION  DATA  AND  COMPUTE  D /DZ 
IF ( ISYM. EQ. 1 )  I SET = 1 

CALL  SMOOTH ( Z i YB 1LZ1JFZ1 JMAX i liNZTLi IU2i ZIi0iDYDXl ) 

CALL  SMOOTH (Zi ZBi LZiJFZi JMAX i liNZTLi IU2iZIi0i0. ) 

1  CONTINUE 

WRITE ( 1 1 200  > 

REWIND  IU2 
DO  6  J= 1 i NZTL 

C  READ  SMOOTHED  DATA  FOR  CONSTANT  2 
CALL  READYZ(J) 

C  SMOOTH  ALONG  Z=CONSTANT  AND  COMPUTE  D/DX 

CALL  SMOOTH ( XBi YBi LX i JFX i IMAXi0iNXTLi IU3i XIi 1 i DUMMY ) 

CALL  SMOOTH ( XBi ZBi LX i JFX i  IMAXi0iNXTLi  IU3 1 X  1 1 2 1 DUMMY ) 

CALL  SMOOTH ( XBi DYBDZ i LX i JFX  i  IMAXi 0i NXTLi  IU3i XI i 4i DUMMY) 
CALL  SMOOTH ( XBi DZBDZ 1LX1JFX1  IMAXi0iNXTLi  IU3i X  1 1 5 i DUMMY ) 
CALL  SMOOTH ( XBi D2YBDZ » LX i JFX  »  I MAX i 0i NXTL i  IU3i XI i 6i DUMMY ) 
CALL  SMOOTH ( XBi D2ZBDZ i LX i JFX i IMAXi0iNXTLi IU3i XIi7iDUMMY) 

C  COMPUTE  HliH2iTHETAiKliK2 
CALL  GEOMET(J) 

WRITE ( 1 i 300 )  J 
REWIND  IU2 
6  CONTINUE 

C  COMPUTE  DERIVATIVES  OF  THETA  ON  EACH  SECTION 
REWIND  IU2 


06100  JFT»NZTL/2+l 

06200  LT  ( 1  )  =  1 

06300  JUP=NZTL-1 

06400  LL-1 

06300  DO  8  J=2»  JUP.2 

06600  LL-LL+ 1 

06700  8  LT(LL)»J 

06800  LT(JFT)=NZTL 

06900  DO  9  I-l.NXTL 

07000  DO  7  J=1.NZTL 

07100  7  THETC J)»THETA( I. J) 

07200  CALL  SMOOTH ( Z I  *  THET »  LT.  JFT » NZTL. 0.  NZTL > IU2. Z I «  8i DUMMY ) 

07300  CALL  SMOOTH ( Z  I .  DTHDZ »  LT »  JFT »  NZTL. 0?  NZTL. IU2. Zlifli DUMMY ) 

07400  9  CONTINUE 

07300  WRITE (1,400) 

07600  REWIND  IU2 

07700  REWIND  IU3 

07800  JFT=NXTL/2+l 

07900  LT  (  1  >  =  1 

08000  IUP=NXTL-1 

08100  LL=1 

08200  DO  13  1=2. IUP.2 

08300  LL=LL+1 

08400  13  LTCLL)=I 

08500  LT ( JFT ) =NXTL 

08600  DO  12  J=1.NZTL 

08700  C  READ  D  (  THETA  > /DZ  FROM  DISK  FILE  FOR  CONSTANT  Z 
08800  CALL  READTH ( J ) 

08900  C  SMOOTH  IHSIA  ALONG  /^CONSTANT  AND  COMPUTE  D/DX 
09000  DO  10  I=1.NXTL 

09100  10  READ(  IU3)  YB<  I ) »  ZB<  I)«H1<I>»  H2<  I )  t  THET  (  I  >  i  CK1  ( I  >  *  CK2  <  I  >. 

09200  CALL  SMOOTH  CXI* THET. LT» JFT. NXTL» 0» NXTL» IU3. XI . 3. DUMMY) 

09300  C  COMPUTE  K12  AND  K21  THEN  WRITE  FINAL  DISK  FILE 
09400  WRITEC IU5. 1200)  ZICJ) 

09300  PSI ( 1 >=PSIZ+H1 C 1 )*DX< 1 > 

09600  DO  11  1=1 . NXTL 

09700  IFC I.LT.NXTL)  PSI ( 1+1 >=PSI ( I )+. 5* (HI ( I )+Hl < 1+1 > ) * < XI < 1+1 )-XI ( I ) ) 

09800  FACT1=CK1 ( I )+DTHDX( I >/Hl ( I ) 

09900  FACT2=CK2 ( I >  +DTHDZ ( I > /H2 ( I ) 

10000  SINTH=S INC  THETA < I, J> ) 

10100  COSTH- COS ( THETA ( I. J) ) 

10200  CK12 ( I >=(— FACTl+COSTH*FACT2> /SI NTH 

10300  CK21 ( I )=(-FACT2+C0STH*FACTl > /SINTH 

10400  WRITEC IU4)  Hi C I ) ,H2< I ) . CK1 ( I > . CK2< I > . CK12C I > . 

10300  *  CK21 ( I ) . THETA C I . J) . PSI C I ) 

10600  WRITEC IU5. 100)  I, XI C I ) , Ye( I ) . ZBC I ) .HI C I ) .H2< I ) , THETA< I . J ) . 

10700  *  CK1 C I ) . CK2C I > . CK12C I ) . CK21 ( I > . PSI < I ) 

10800  11  CONTINUE 

10900  REWIND  IU2 

11000  12  CONTINUE 

11100  100  FORMAT ( I4.3X. 1P11E11.3) 

11200  200  FORMAT ( 1 X .’ A 1 1  Section  Data  Smoothed  &  Differentiated’//) 

11300  300  FORMAT  C  1  X .’ Geometr  i  ca  1  Parameters  Computed  for  J  «=  ’.13) 

11400  400  FORMAT  C / 1 X. ’ d  C  Theta ) /dz  Computed  for  all  Sections’/) 

11500  1000  FORMAT C3F10. 5. 11X.2F10.3) 

11600  1100  FORMATC 10X.2F10.S. 11X.2F10. 5) 

11700  1200  FORMAT C // 1 X » ’ Z  =  ’ *F7. 3/3X1HI . 10X2HXB. 9X2HYB. 9X2HZB, 

11800  *  9X2HH1 . 9X2HH2. 8X5HTHETA. 7X2HK1 . 9X2HK2. 

11900  *  9X3HK12. 8X3HK21 . 0X3HPSI ) 

12000  END 
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SUBROUTINE  ANGLE  <  COSTH*  SINTH, THETA) 

******************************** ************************************ 
THIS  SUBROUTINE  CALCULATES  AN  ANGLE  IN  EITHER  THE  FIRST  OR  SECOND 
QUADRANT  GIVEN  THE  COSINE  OF  THE  ANGLE 
******************************************************************** 

IF < COSTH . EQ- 0 . )  GO  TO  2 
IF  <  COSTH. GT. 0.  )  GO  TO  1 
SINPHI =— COSTH 
COSPHI=SQRT ( 1 . -SINPHI**2> 

TANPH I =S I NPH I / COSPH I 
PH I =AT AN ( T AN PH I ) 

THETA=355 • /226. +  PH I 
SI NTH=COSPH I 
RETURN 

SINTH-SQRT< 1 . -C0STH**2> 

TANTH=S I NTH / COSTH 
THETA=AT AN <  TANTH ) 

RETURN 
S I NTH= 1 . 

THETA=355 • /226. 

RETURN 
END 

SUBROUTINE  ARC ( DYDX 1 ) 

********************************************************  ************-. 
THIS  SUBROUTINE  COMPUTES  LOCAL  ARCLENGTH  AND  RESCALES  TO  UNITY 
********************************************************************.; 
COMMON/ COORD/  Z ( 40 ) 

COMMON/ COORDB/XB ( 40 ) *  YB<40>  *  ZB (40) 

COMMON/NDATA/  IUl ,  IU2,  IU3*  I U*  *  IU5,  IMAX, JMAX • NYTL« NZTL, J«=*Z*LZ (30) 
Z<1>^0. 

DO  1  J=2, JMAX 
JM=J-1 

DYB=YB  <  J ) — YB <  JM ) 

DZB=ZB ( J ) — ZB  <  JM ) 

DZ-SQRT ( DYB*DYB+DZB*DZB ) 

Z ( J)=Z ( JM)+DZ 
ZSCALE=Z( JMAX ) 

DO  2  Jsa2 1  JMAX 
Z( J)=Z( J)/ZSCALE 
DYDX 1 ZSCALE 
RETURN 

END  # 

SUBROUTINE  GEOMET(J) 

********************************************************************** 
THIS  SUBROUTINE  CALCULATES  THE  METRIC  AND  CURVATURE  PARAMETERS 
****************  **************************************************** -» 
COMMON/ COO RDB/ XB  <  40 ) ?  YB ( 40 ) i ZB( 40 ' 

COMMON/GEOMP/  HI <  40 ) *  H2 ( 40 ) » THETA (40, 40) , CK1 (40) , CK2(40) ,  • 

*  CK12(40) t CK21 (40) , PSI (40) 

COMMON/TEMPZ /  DYBDZ ( 40 ) , D2YBDZ (40) , DZBDZ (40 ) , D2ZBDZ ( 40) 

COMMON/ TEMP X /  DYBDX (  40 )  , D2YBDX (40), DZBDX (40), D2ZBDX ( 40 ) 

COMMON/NDATA/  IUl* IU2, IU3, IU4, IU5, I MAX , JMAX , NXTL , NZTL i JFZ , LZ ( 30 ) 

DO  3  1=1, NXTL 

HI < I )=SQRT< 1 . +DZBDX ( I )**2+DYBDX ( I )**2)  • 

H2< I )=SORT( DYBDZ ( I )**2+DZBDZ( I )**2) 

COSTH=( DYBDZ ( I ) *DYBDX ( I )+DZBDX ( I )*DZBDZ ( I ) ) / (HI ( I ) *H2 ( I ) ) 

CALL  ANGLE ( COSTH, SI NTH, THETA ( I, J) > 

FACT 1=DYBDZ ( I )*DZBDX( I ) -DYBDX ( I )*DZBDZ( I) 

FACT2-DZBDX ( I ) *D2YBDX ( I ) -DYBDX ( I ) *D2ZBDX ( I ) 

FACT3=DYBDZ ( I ) «D2YBDX ( I ) +DZBDZ < I > *D2ZBDX ( I )  • 
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18100  DEN0M=H1 ( I )**4*H2< I )*SINTH 

18200  CK 1 < I ) - ( FACT 1 *FACT2+FACT3 ) /DENOM 

18300  FACT2=DZBDZ( I )*D2YBDZ< I ) -DYBDZ ( I >*D2ZBDZ ( I ) 

18400  DEN0M=H1 ( I )*H2( I )**4*SINTH 

18300  CK2< I >=-FACTl*FACT2/DEN0M 

18600  WRITE< IU3)  YB< I > i ZB< I ) ,H1 ( I > » H2< I ) , THETA < I . J ) » CK1 < I ) i CK2^ I ) 

18700  3  CONTINUE 

18800  RETURN 

18900  END 

19000  SUBROUTINE  GRID 

19100  C  ********************** ********* **************************** ********* 
19200  C  THIS  SUBROUTINE  SETS  UP  THE  COMPUTATIONAL  GRID 

19400  COMMON/ INTRP/  X  I ( 40 ) » Z I < 40) 

19500  COMMON/ X DATA/  PS I Z . XSTART i XEND , DX ( 1 0 ) , NX < 1 0 > 

19600  COMMON/ZDATA/  Z START •>  ZEND i  DZ  (  10)  »NZ  (  10) 

19700  COMMON/NDATA/  IU1 , IU2i IU3. IU4, IU5 i I MAX » JMAX » NXTL i NZTL » JFZ » LZ < 30 ) 

19800  X  I ( 1 ) =XSTART+DX ( 1 ) 

19900  K=1 

20000  MCUT =NX ( 1 ) 

20100  NXTLM=NXTL— 1 

20200  DO  1  1  =  1 1  NXTL.M 

20300  IF( I . LE. MCUT)  GO  TO  1 

20400  K=K+ 1 

20500  MCUT =MCUT +NX  <  K ) 

20600  1  XI < 1+1 )*XI ( I )+DX(K) 

20700  IF ( ABS ( X  I <  NXTL ) —XEND ) . GT.  .  1*DX ( K ) )  GO  TO  7 

20800  XI (NXTL)=XEND 

20900  Z I < 1 ) =ZSTAPT 

21000  K=1 

21100  MCUT =NZ  < 1 ) 

21200  NZTLM=NZTL- 1 

21300  DO  3  1  =  1*  NZTLM 

21400  IF( I.LE.MCUT)  GO  TO  3 

21500  K=K+1 

21600  MCUT=MCUT+NZ ( K ) 

21700  3  ZI ( 1+1 >=ZI ( I )+DZ (K) 

21800  IF ( ABS ( Z I ( NZTL )  —  ZEND ) . GT. .  1*DZ<K> )  GO  TO  8 

21900  Z I ( NZTL ) =ZEND 

22000  DO  4  1=1* NXTL 

22100  4  WR I TE  (  I U4  )  XKI) 

22200  DO  5  1=1* NZTL 

22300  5  WRITE ( IU4  >  ZI(I) 

22400  RETURN 

22500  7  WRITE< IU5,  100)  X  I ( NXTL ) * XEND 

22600  STOP 

22700  8  WRI TE ( I U5  *  200  >  Z I ( NZTL > . ZEND 

22800  STOP 

22900  100  FORMAT( /IX* ’STEPSIZE  ARRAY  IMPLIES  XEND  =’*1PE13.6* 

23000  *’  WHICH  DOES  NOT  MATCH  THE  DESIRED  VALUE _ ’*E13.6) 

23100  200  F0RMAT(/1X*  ’STEPSIZE  ARRAY  IMPLIES  ZEND  *  1 PE13. 6* 

23200  *’  WHICH  DOES  NOT  MATCH  THE  DESIRED  VALUE. . . ’ * E13. 6 ) 

23300  END 

23400  SUBROUTINE  IN I TL ( X  * Y , N* L * JF,  I START * DYDX 1 ) 

23500  C  ****** ***#*-**■****•*•*■*-* *■*#**•**■*•* ************************************** 
23600  C  THIS  SUBROUTINE  SETS  UP  INITIAL  VALUES  FOR  THE  VARIOUS  ARRAYS 
23700  C  ******************************************************************** 
23800  DIMENSION  X ( 40 )iY( 40 )iL( 30 ) 

23900  COMMON/OUAN/  A ( 30 ) i B ( 30 ) i PH0PH0 < 30 ) , PH0PH 1 < 30 ) , 

24000  *  PH0PS0 ( 30 ) i PH0PS 1 <  30 ) »  PHI  PHI <  30 ) »  PHI PS0 ( 30 ) » 


24100 
24200 
24300 
24400 
24500 
24600 
24700 
24800 
24900 
25000 
25100 
25200 
25300 
25400 
25500 
25600 
25700  2 

25800 
25900 
26000 
26100 
26200 
26300 
26400 
26500 
26600 
26700  3 

26800 
26900 
27000 
27100 
27200 
27300 
27400 
27500 
27600 
27700 
27800 
27900 
28000 
28100 
28200 
28300 
28400 
28500 
28600 
28700 
28800  1 
28900 
29000 
29100 
29200  C 
29300  C 
29400  C 
29500  C 
29600 
29700 
29800 
29900 
30000 


*  PH1PS1 (30) »  PI (30) , P2<30>  *  P3<30) »Q1 (30) »  ~  V] 

*  Q2 ( 30 ) *  03 ( 30 ) i X  I ( 30 ) i YPHI0<30> i YPHI 1 ( 30 ) t  j 

*  H ( 30 ) 

COMMON/DIMN/  JFM 
COMMON/LDATA/  ISYfit  ISET 
JFM=JF-1 
Y 1  =Y  (  1  ) 

YN«Y(N) 

IF ( I START, NE. 0 )  GO  TO  2 
X32— X ( 3 ) “X ( 2 ) 

X31  =  X ( 3 ) — X  < 1 ) 

X21 =X  <  2 ) —X ( 1 ) 

AA=-<X21+X31 >/<X21*X31 ) 

BB“X31 / < X21*X32> 

CC=-X21/<X31*X32> 

DYDX  1  =AA*Y  (  1  >+BB*Y<2)+CC*Y<3> 

X32=X ( N—2 ) —X ( N— 1 ) 

X31«X(N-2>-X(N> 

X2 1  =  X  <  N— 1 ) —X  <  N ) 

AA==—  (  X21  +  X31  >/<X21*X31  > 

BB—X31 / <  X21*X32 ) 

CC=-X21/<  X31*X32) 

DYDXN=AA*Y  <  N  >  -+  BB*Y  (  N- 1  )  +  CC*Y  (  N-*2  ) 

IF ( ISET- EQ. 0 )  GO  TO  3 
DYDXN=0. 

YN=~ <  BB*Y <  N- 1  ) +CC*Y <  N-2 )  ) /AA 
ISET =0 
PI ( 1 )«Y1 
P2 ( i >  =0. 

P3< 1 >«0. 

0.1  <  1  )  =  DYDX  1 

02.  (  1  >=0. 

03 < 1 ) =0 • 

A  <  JF ) =YN 
B ( JF ) =DYDXN 
DO  1  J= 1 »  JF 
K=L ( J ) 

XI < J)=X<K) 

PH0PH0 < J)=0. 

PH0PH1 < J)=0. 

PH0PS0  <  J ) =0 - 
PH0PS1 < J)^0. 

PHI  PHI ( J)=0. 

PHI PS0  <  J ) “0- 
PHI  PS 1 < J)=0. 

YPHI 0  <J>*=0, 

YPHI1  <J)=0. 

CONTINUE 
RETURN 
END 

SUBROUTINE  I NTERP ( I i J i XXi IU2» I FLAG ) 

********************** ***<*»**^ ************ ******************* ****•»+ 


THIS  SUBROUTINE  USES  THE  SPLINE  FIT  TO  INTERPOLATE  THE  ORIGINAL 

DATA  AND  TO  COMPUTE  ITS  FIRST  AND  SECOND  DERIVATIVES  • 

s 

********************************************************************* 
COMMON/QUAN/  A ( 30 )i6( 30 ) » PH0PH0 ( 30 ) , PH0PH 1 ( 30 ) » 

*  PH0PS0 ( 30 ) i PH0PS1 (30) »  PHI  PHI ( 30 ) , PH1PS0(30) * 

*  PHI  PS  1  <  30  )  »  P 1  (  30  )iP2(  30  )iP3(  30  )  i  <3 1  (  30  )  * 

*  Q2(30) i Q3<30> i XI (30) . YPHI0C30) » YPHI 1 (30)  . 

*  H ( 30  >  • 
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30100 

30200 

30300 

30400 

30500 

30600 

30700 

30000  1 

30900 

31000 

31100 

31200 

31300 

31400 

31500 

31600 

31700 

31800 

31900 

32000 

32100 

32200 

32300 

32400 

32500 

32600 

32700 

32800 

32900  2 

33000 

33100 

33200 

33300 

33400 

33500  3 

33600 

33700 

33800 

33900  4 

34000 

34100 

34200 

34300 

34400 

34500 

34600 

34700 

34800  C 

34900  C 

35000  C 

35100 

35200 

35300 

35400 

35500 

35600 

35700 

35600 

35900 

36000  1 


COMMON/ TEMP X /  DYBDX ( 40  >  t  D2YBDX (40)i DZBDX <40)» D2ZBDX C  40 ) 
COMMON/TEMPZ /  DYBDZ (40)i D2YBDZ (40)i DZBDZ (40)» D2ZBDZ ( 40 ) 

COMMON/ COORDB/XB( 40 ) t  YB  ( 40) t ZB (40) 

COMMON/DTHET/  DTHDX <  40 ) t  DTHDZ ( 40 )  t THET ( 40  >  t  LT ( 30  > 

COMMON/D IMN/  JFM 
XFCXX.LE. XI < 1+1 > >  GO  TO  1 
I®I  +  1 
IP=I+1 

DX« ( XX  — X  I ( I ) >/H< I ) 

DXP= (XI ( IP)-XX  >/H( I ) 

PHI 1  =  (3. -2. *DX )*DX**2 
PHI0=1 . —PHI 1 
PSI 1--H( I )*DXP*DX**2 
PSI0=H< I )*DX*DXP**2 
PHI  lP=-6.  *DX*(  1  •  — DX  )  /H (  I  ) 

PHI0P=— PHI 1 P 

PSI 1 P=-DX* ( 2. *DXP-DX ) 

PSI0P=DXP*(DXP-2. *DX ) 

PH1PP=6. *( 1.-2. *DX) /H( I )**2 
PH0PP— -PHI PP 

PS1PP=2. *(2. *DX— DXP ) /H ( I ) 

PS0PP=2.*(DX-2.*DXP)/H( I ) 

YY=A ( I ) *PH I 0  +A(IP)*PHI1  +B(I)*P5I0  +B(IP)*PSI1 
YP=A  (  I  )  *PHI0P+A  (  I  P  )  *PHI  1  P+B  (  I  )  *PS I0P+B  (  IP)*PSI1P 
Y2=A ( I ) *PH0PP+A ( I P ) * PH 1 P P+B  < I ) *PS0PP+B( I P ) *PS1 PP 
IF ( I FLAG. GT. 0 )  GO  TO  2 
WRIT£( IU2)  YYt YPt Y2 
RETURN 

ih l 1KLAG. GT. 2)  GO  TO  4 
I F ( I FLAG . EQ . 2 )  GO  TO  3 
YB  ( J )  =*YY 

DYBDX (J)^YP 
D2YBDX ( J ) =Y2 
RETURN 
ZB ( J ) ~YY 
DZBDX ( J ) =YP 
D2ZBDX ( J ) ~Y2 
RETURN 

IF ( I FLAG. EQ.  3)  DTHDX (J)=YP 
IF ( I FLAG. EQ. 4 )  DYBDZ ( J ) =YY 
IF ( I FLAG. EQ. 5 )  DZBDZ (J)=YY 
I F  < I FLAG • EQ . 6 )  D2YBDZ(J)=YY 
IF ( I FLAG. EQ. 7)  D2ZBDZ <  J ) =YY 
IF ( I FLAG. EQ. 8 )  DTHDZ(J)=YP 
RETURN 
END 

SUBROUTINE  NAM IN 

********************************************************************* 
THIS  SUBROUTINE  COORDINATES  READING  OF  INPUT  DATA 
********************************************************************* 
COMMON/NDATA/  N ( 40 ) 

COMMON/MDATA/  M(31) 

COMMON/LDATA/  ISYMt ISET 
COMMON /X DATA/  X(13)tNX(10) 

COMMON/ZDATA/  Z( 12) tNZ( 10) 

ISYM=0 

CALL  NAMLST ( 6t  2t  ISYMt  DUMMY ) 

DO  1  I»ltl0 
N< I >«0 

CALL  NAMLST (6t2t N( I ) t DUMMY) 
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36100 

36200 

36300 

36400 

36500 

36600 

36700 

36800 

36900 

37000 

37100 

37200 

37300 

37400 

37500 

37600 

37700 

37800 

37900 

38000 

38100 

38200 

38300 

38400  C 

38500  C 

38600  C 

38700 

38800 

30VW0 

39000 

39100 

39200 

39300 

39400 

39500 

39600 

39700 

39800 

39900  C 

40000  C 

40100  C 

40200 

40300 

40400 

40500 

40600 

40700 

40800 

40900 

41000 

41100 

41200 

41300 

41400 

41500 

41600  C 

41700  C 

41800  C 

41900  C 

42000 


3 

4 


5 

6 


7 


1 


2 

3 


IUP=10+N< 10) 

DO  2  1  =  11, IUP 
N< I  )«0 

CALL  NAMLST<6,2,N< I ) ,DUMMY> 

M< 1  )=0 

CALL  NAMLST  <6i2iM( 1 )t DUMMY ) 

IUP=1+M< 1 ) 

DO  3  1=2, IUP 
M ( I ) =0 

CALL  NAMLST  <6»2iM( I ) i DUMMY ) 

DO  4  1  =  1  *  13 

CALL  NAMLST < 6,2, li  X< I ) ) 

DO  5  1=1 i 10 
NX ( I ) =0 

CALL  NAMLST  <6*2*  NX  < I ) *  DUMMY ) 

DO  6  1=1, 12 

CALL  NAMLST (6,2,  1 ,  Z  (  I  )  ) 

DO  7  1=1, 10 
NZ ( I ) =0 

CALL  NAMLST (6,2, NZ ( I ) , DUMMY ) 

RETURN 
END 

SUBROUTINE  NAMLST ( 1 1 ,  I O, N, X ) 
******************************************************************** v 
THIS  SUBROUTINE  READS  INPUT  DATA 


********************************************************************- 
DIMENSION  A  <  3 ) 

IF ( N  .NE.  0)  GO  TO  1 

A  d  >  ,  A  C  2  >  ,  A  <  3  >  ,  N 
A  < 1 > , A<2) , A<3> , N 


A< 1 ) , A  <  2 ) i A ( 3 ) ♦ X 
Ad  >  f  A<2>!  A<3>  «  X 


KtiAu  111*2) 

WRITE < 10,2) 

RETURN 
READ  <11, 3) 

WR I TE  < 10*3) 

RETURN 

FORMAT  < IX, 3A4 ,  14) 

FORMAT  < IX, 3A4 , 1  PE 13. 6) 

END 

SUBROUTINE  READTH(J) 

********************************************************************- 


THIS  SUBROUTINE  READS  D< THETA )/DZ  ON  A  SECTION  FROM  DISK 
********************************** *********************************** 
COMMON/DTHET/  DTHDX <  40 ) , DTHDZ <  40 ) , THET  <  40 ) , LT  <  30 ) 

COMMON/NDATA/  IU1, IU2, IU3, IU4, IU5, I MAX , JMAX , NXTL , NZTL , JFZ , LZ < 30 ) 

DO  3  1=1, NXTL 
DO  1  JJ=1,J 

1  READ  < IU2 )  YY , D, D 

DTHDZ < I >=YY 
IF < J. EQ. NZTL )  GO  TO  3 
JP= J+l 

DO  2  JJ=JP,NZTL 

2  READ  < IU2 )  D,D,D 

3  CONTINUE 
RETURN 
END 

SUBROUTINE  READYZ  <  J ) 
***************************************************  ****************** 
THIS  SUBROUTINE  READS  SMOOTHED  SECTION  DATA  FROM  DISK  AND  FILLS 
ARRAYS  IN  PREPARATION  FOR  SMOOTHING  IN  THE  STREAMWISE  DIRECTION 
*************************************  *******************************  * 
COMMON/ COORDB/XB  <  40  > , YB<40> , ZB<40) 


42100 

42200 

42300 

42400 

42500 

42600 

42700 

42800 

42900 

43000 

43100 

43200 

43300 

43400 

43500 

43600 

43700 

43800 

43900 

44000 

44100 

44200 

44300 

44400 

44500 

44600  C 

44700  C 

44800  C 

44900  C 

45000 

45100 

45200 

45300 

45400 

45500 

45600 

45700  C 

45800 

45900 

46000  C 

46100 

46200 

46300  C 

46400 

46500 

46600 

46700 

46800 

46900 

47000 

47100 

47200 

47300 

47400 

47500 

47600 

47700 

47800  C 

47900 

48000 


COMMON/TEMPZ /  DYBDZ <  40  > »  D2YBDZ  <  40 ) *  DZBDZ ( 40 ) * D2ZBDZ ( 40 ) 
COMMON/NDATA/  IU1  * IU2* IU3* IU4* IU5* I MAX i JMAX » NXTL * NZTL * JFZ * LZ < 30 ) 

DO  6  1=1* IMAX 
DO  1  JJ= liJ 

1  READ ( IU2 )  YYiYP*Y2 
YB< I)=YY 

DYBDZ ( I >=YP 
D2YBDZ ( I >=Y2 
IF< J.EQ.NZTL)  GO  TO  3 
JP=J+ 1 

DO  2  JJ=JP*NZTL 

2  READ  < IU2 )  D*D*D 

3  CONTINUE 

DO  4  JJ=1*J 

4  READ ( IU2 )  ZZ *  ZP*  Z2 
ZB < I )=ZZ 

DZBDZ ( I ) =ZP 
D2ZBDZ  < I  )=Z2 
IF< J.EQ.NZTL)  GO  TO  6 
DO  5  JJ= JP*  NZTL 

5  READ ( I U2 )  D*D*D 

6  CONTINUE 
RETURN 
END 

SUBROUTINE  SMOOTH < X * Y* L * JF* N * ISTARTiNN,  IU2*  XT*  IFLAG*  DYDX 1 ) 
******************************************  *************************** 
THIS  SUBROUTINE  FITS  A  TABULATED  FUNCTION  WITH  CUBIC  SPLINES  USING 
LEAST-SQUARES  ERROR  MINIMIZATION 

•a-**************************************************************-***-**..* 

DIMENSION  X ( 40 ) *  Y ( 40 ) *  L  < 30 ) »  XT ( 40 ) 

COMMON/QUAN/  A  <  30 ) *  B  <  30 ) *  PH0PH0 ( 30 ) , PH0PH1 ( 30 ) * 

*  PH0PS0 ( 30 ) * PH0PS1 (30) *  PHI  PHI (30) *  PH1PS0(30) * 

*  PH1PS1 (30) *  PI (30) *  P2(30) *  P3(30) *  Q1 (30) * 

*  02(30) *03(30) *  XI (30) * YPHI0(30) *  YPHI1 ( 30 ) * 

*  H ( 30 ) 

COMMON/ D I MN/  JFM 

INITIATE  THE  COMPUTATION 

CALL  INI TL ( X  *  Y *  N *  L *  JF *  I START* DYDX 1 ) 

CALL  SUMUP ( X  *  Y*  N*  L  *  JF ) 

COMPUTE  THE  P’S  AND  Q’S 
DO  1  J=2  *  JFM 
M=J—  1 

COEFFICIENTS  IN  THE  ORIGINAL  EQUATIONS 
A 1 =PH0PH 1 ( M ) 

B1=PH1PH1 ( M ) +PH0PH0  <  J ) 

C 1 =PH0PH 1 < J) 

D1=PH1PS0(M) 

E1=PH1PSI (M ) +PH0PS0 ( J ) 

F1«PH0PS1 < J) 

G1=YPHI 1  ( M ) +YPH I 0  (J) 

A2=3. /H(M)**2 

B2=3. *  < 1 . /H ( J ) **2—1 • /H  <  M ) **2 ) 

C2=— 3. /H< J>**2 
D2=l. /H(M) 

E2=2.#(  1 .  /Hi  J>  +  1.  /H(M>  ) 

F2=l. /H(J) 

G2*=0. 

MODIFIED  COEFFICIENTS 

B1S=B1+A1 *P2  <  M  >+D 1*02 (M  > 

E1S=E1+A1*P3(M)+D1*Q3(M) 


48100 

48200 

48300 

48400 

48300  C 

48600 

48700 

48800 

48900 

49000 

49100 

49200 

49300  i 

49400  C 

49500 

49600 

49700 

49800 

49900 

50000 

50100 

50200  2 

50300  C 

50400 

50500 

50600 

50700  2 

50800 

-J  HJ  / 

51000 

51100  C 

51200  C 

51300  C 

51400  C 

51500 

51600 

51700 

51800 

51900 

52000 

52100 

52200 

52300 

52400 

52500 

52600 

52700 

52800 

52900 

53000 

53100 

33200 

53300 

33400 

53300 

33600 

53700 

53800 

33900 

34000 


G1S=G1-A1*P1 (M)-D1*01 <M> 

B2S=B2+A2*P2  <  M ) +D2+02 ( M ) 

E2S=E2+A2*P3  <  M )  +D2*G3  ( M ) 

G2S=G2-A2*P1 (M)-D2*01 <M> 

THE  P’S  AND  O’S 

DEL= 1 . / ( B 1 S*E2S-B2S * E 1 S ) 

PI ( J>=(E2S*G1S-E1S*G2S)*DEL 
P2< J)=(E1S*C2  -E2S*C1  >*DEL 
P3< J)=(E1S*F2  — E2S*F 1  )*DEL 
01 < J)=(B1S*G2S-B2S*G1S)*DEL 
G2( J)=(B2S*C1  -B1S*C2  )*DEL 
Q3( J)=(B2S*F1  -B1S*F2  >*DEL 
CONTINUE 

SOLVE  FOR  THE  A’S  AND  B’S 
K=JFM 
KP=JF 

DO  2  J=1 »  JFM 

A  <  K ) -P 1  ( K )  +  P2 ( K ) *A  <  KP ) +  P3 ( K ) *B  < KP ) 

B  <  K  >  =0 1 ( K ) +02 ( K ) * A  <  KP ) +Q3  <  K ) *B ( KP ) 

KP=K 
K=K- 1 
CONTINUE 

INTERPOLATE  TO  DESIRED  MESH 
1  =  1 

DO  3  J= 1 i NN 

CALL  INTERP( I, J. XT< J) , IU2» IFLAG) 

CONTINUE 

RETURN 

END 

SUBROUTINE  SUMUP ( X » Yi N» L . JF ) 

*******»*************4Ht*************»**X«-*********************.|HH{..»*« 

THIS  SUBROUTINE  PERFORMS  SUMMING  OPERATIONS  ON  THE  CUBIC-SPLINE 
WEIGHTING  FUNCTIONS 

DIMENSION  X ( 40 ) i Y ( 40 ) »  L ( 30 ) 

COMMON/OUAN/  A ( 30 >  »  B  (  30 ) » PH0PH0 ( 30 )  »  PH0PH1 ( 30 ) » 

*  PH0PS0 ( 30 ) »  PH0PS1 (30) i PHI  PHI (30) »  PH1PS0(30) t 

*  PH1PS1 ( 30 ) i PI ( 30 ) »  P2 ( 30 ) i P3 ( 30 ) iQl ( 30 ) < 

*  02(30) »  03(30) >  XI (30) i YPH 10 ( 30 ) t YPHI 1 ( 30 > i 

*  H ( 30 ) 

COMMON/DIMN/  JFM 
DO  2  J= 1 i JFM 

H( J ) =X I ( J+l ) —X  I ( J  > 

KL»L( J) 

KU=L( J+i ) 

DO  1  K=KL »  KU 
DX=(X(K>— XI (J) )/H(J) 

DXP=( XI ( J+l >-X(K> ) /H ( J ) 

YY=Y  <  K ) 

PHI 1=(3. -2. *DX )*DX««2 
PH 1 0=1 . -PHI  1 
PS I 1«-H( J)*DXP*DX**2 
PSI0=  H( J)*DX#DXP+*2 
PH0PH0 ( J  >  =  PH0PH0 ( J ) +  PH 1 0*  PH I  0 
PH0PH1 ( J)=PH0PH1 ( J)+PHI0*PHI 1 
PH0PS0 ( J ) =PH0PS0 ( J ) +PHI 0+PS I 0 
PH0PS 1 ( J ) =PH0PS 1 ( J)+PHI0*PSI 1 
PHI  PHI ( J)=PH1PH1 ( J)+PHI 1*PHI 1 
PH 1 PS0 ( J ) =PH 1 PS0 ( J ) +PH I  1  *  PS 1 0 
PHI  PS 1 ( J  >«PH1PS1 ( J)+PHI l  +  PSI 1 
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00100  PROGRAM  VELOC 

00200  C  **#*#***•*#**#*****#****##*******#****##****##**********************• 
00300  C  THIS  PROGRAM  SMOOTHS  AND  INTERPOLATES  BOUNDARY  LAYER  EDGE  VELOCITY 
00400  C  COMPONENTS  ONTO  A  RECTANGULAR  GRID  IN  TRANSFORMED  COORDINATES  AND 
00500  C  WRITES  A  DISK  FILE  FOR  USE  BY  PROGRAM  EDDY3 

00600  C  *#**♦*##** w*******-*-*************-*****-*-* ************* **************** 
00700  COMMON/COORD/  Z<40> 

00800  COMMON/ COORDB/ XB ( 40  > i YB ( 40 ) i ZB ( 40 ) 

00900  COMMON/DERI V/  DYBDZ ( 40 ) i DZBDZ ( 40 ) i DWDZ 1 ( 40 ) *  DWDZN ( 40 ) 

01000  COMMON/ INTRP/  X  I < 40 ) , Z I ( 40 ) 

01100  COMMON/VEL/  UE < 40 > » VE < 40 ) i WE < 40 ) 

01200  COMMON/VEL I /  UEI ( 40. 40 ) t WEI ( 40i 40 ) 

01300  COMMON/NDATA/  IU1 * IU2. IU3» IU4 » IU5 » IMAX » JMAX » NXTL » NZTL > JFZ » LZ < 30 ) 

01400  COMMON/MDATA/  JFX»LX<30> 

01500  COMMON/LDATA/  ISYM» ISET 

01600  CALL  OPEN  ( 6»  ’  NAMVEl.  M) 

01700  CALL  NAM IN 

01800  CALL  OPEN < IU1 1  ’ UEWE  ’,1> 

01900  CALL  OPEN < IU3» ’VELOCITY  ’>8) 

02000  CALL  OPEN ( I U4, ’ GEOMETRY  ’.32) 

02100  DO  1  1=1 , NXTL 

02200  1  READ ( IU4  >  XI (I) 

02300  DO  2  1  =  1  *  NZTL 

02400  2  READ ( IU4 )  ZI(I) 

02500  JMAX=JMAX+2 

02600  JMAXM=JMAX— 1 

02700  DO  4  1=1 i IMAX 

02800  ISET =0 

S>2Y00  C  HEAD  St  C  I  I  ON  VELOCITIES 

03000  READ ( I U 1 1 1 000 )  XB ( I ) »  YB (2) i ZB(2) i CX »  CY »  CZ  »  Q 

03100  UE  ( 2 )  =<5*CX 

03200  VE  <  2 ) =0*CY 

03300  WE  <  2  >  =0*CZ 

03400  DO  3  J=3»  JMAXM 

03500  READ( IU1 » 1 100)  YB( J) » ZB( J) » CXi CYt CZ.O 

03600  UE(J)=0*CX 

03700  VE<J)=(5*CY 

03800  3  WE(J)*Q*CZ 

03900  C  COMPUTE  VELOCITIES  ON  KEEL  AND  FREE  SURFACE 
04000  CALL  ENDS 

04100  C  COMPUTE  AND  RESCALE  LOCAL  ARCLENGTH 
04200  CALL  ARC ( DYDX 1 ) 

04300  C  COMPUTE  TANGENTIAL  VELOCITY  ON  SECTION 
04400  IF( ISYM.EO. 1 )  ISET=1 

04500  CALL  SMOOTH ( Z  *  YB iLZi JFZ »  JMAX i  li JMAX 1Z161 DYDX 1 1 I DUM ) 

04600  CALL  SMOOTH  <  Z i ZBi LZ i JFZ  »  JMAX » 1 1 JMAX i Zi7i0. i  I DUM ) 

04700  DO  10  J=1 »  JMAX 

04800  H2-S0RT  <  DYBDZ  <  J ) **2+DZBDZ ( J  >  **2 ) 

04900  Q2=UE  <  J  )  **2+' VE  (  J  )  **2+WE  <  J  >  **2 

05000  WE< J)=(VE( J)*DYBDZ( J ) +WE ( J>*DZBDZ< J) ) /H2 

05100  10  UE  <  J ) "SORT ( 02— WE ( J ) **2 ) 

05200  C  SMOOTH  AND  INTERPOLATE  ORTHOGONAL  VELOCITIES 
05300  CALL  SMOOTH ( Z  »  UE» LZ  »  JFZ  »  JMAX  *  l»NZTLtZIi0»0. t I) 

05400  CALL  SMOOTH ( Z  »  WE iLZi JFZ i JMAX  »  0  f  NZTL  t  Z 1 1 1 «  DUMMY » I ) 

05500  WRITE ( 1 i 150 )  I 

05600  4  CONTINUE 

05700  WRITE (It  200 ) 

05800  C  SMOOTH  AND  INTERPOLATE  DATA  ALONG  Z=CONSTANT 
05900  DO  6  J-l.NZTL 

06000  DO  5  1*1  * IMAX 
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06100 

06200 

06300 

06400 

06500 

06600 

06700 

06800 

06900 

07000 

07100 

07200 

07300 

07400 

07500 

07600 

07700 

07800 

07900 

08000 

08100 

08200 

08300 

08400 

08500 

08600 

08700 

08800 

66900 

09000 

09100 

09200 

09300 

09400 

09500 
09600 
09700 
09800 
09900 
10000 
10100 
10200 
10300 
10400 
10500 
10600 
10700 
10800 
10900 
11000 
11100 
11200 
11300 
11400 
1  1500 
11600 
11700 
11800 
11900 
12000 


UE  < I ) =UE I  <  I ,  J  ) 

5  WE ( I ) »WE I ( 1 9  J ) 

CALL  SMOOTH ( XB,  UEi LX ,  JFX , I MAX , 0  9  NXTL  tX!i2i DUMMY »  IDUM) 

CALL  SMOOTH (XB»WEi LX i JFX, IMAX,0,NXTL, XI, 3, DUMMY, IDUM) 

C  COMPUTE  NONORTHOGONAL  VELOCITY  COMPONENTS 
CALL  ROTATE ( J ) 

6  CONTINUE 

CALL  SMOOTH <XB,DWDZ1, LX, JFX, IMAX,0,NXTL, XI, 4, DUMMY, IDUM) 

CALL  SMOOTH  <  XB , DWDZN* LX,JFX,IMAX,0,NXTL,XI,5, DUMMY , I DUM ) 

C  WRITE  DISK  FILE  FOR  SMOOTHED  DATA  ALONG  Z=CONSTANT 
DO  7  J— 1 , NZTL 
WRITE ( IU5 , 1200)  Z I < J ) 

DO  8  1=1, NXTL 

WRITE ( I U3 )  UEI ( I, J) , WEI <  I,o  ) 

8  WRITE< IU5 ,  100)  I, XI < I ) ,UEI < I , J ) , WE  X ( 1 1 J  > 

7  CONTINUE 

WR I TE ( IU5, 250) 

DO  20  1=1, NXTL 

WRITE ( IU3 )  DWDZ1 < I ) , DWDZN ( I ) 

20  WRITE( IU5,  100)  I , X  I < I ) , DWDZ 1(1), DWDZN ( I ) 

100  FORMAT ( 14, 5X, 1P3E11.3) 

150  FORMAT ( IX, ’Section  Data  Smoothed  &  Differentiated  for  I  =’,I3) 

200  FORMATC IX,  1  All  Section  Data  Smoothed  &  Differentiated’//) 

250  FORMAT  < / /3X 1HI , 10X2HXB, 8X5HDWDZ 1 , 7X5HDWDZN ) 

1000  FORMAT ( 3F1 0. 5 , 4F 10.6) 

1100  FORMAT ( 1 0X , 2F 10.5, 4F 10.6) 

1200  FORMAT (//IX,* Z  =  9 , F7. 3/3X1HI, 10X2HXB, 9X2HUE, 9X2HWE ) 

END 

SuBkOUTInl  akCU)yl)xi  ) 

C  ************************************************************************* 
C  THIS  SUBROUTINE  COMPUTES  LOCAL  ARCLENGTH  AND  RESCALES  TO  UNITY 
C  *******************•*■****************************-**#**#****■****#**•*.#.** 
COMMON/ COORD/  Z < 40 ) 

COMMON/ COO RDB/ XB ( 40 ) , YB ( 40 ) , ZB (40) 

COMMbN/NDATA/  IU1 , IU2, IU3, IU4, IU5, I MAX , JMAX , NXTL, NZTL, JFZ , LZ ( 30 ) 

Z( 1 )=0- 

DO  1  J=2,JMAX 

JM= J— 1 

DYB=YB< J)-YB( JM) 

DZB=Z3< J)-ZB< JM) 

DZ=SORT ( DYB*DYB+DZB*DZe ) 

1  Z< J)=Z< JM)+DZ 
ZSCALE=Z< JMAX ) 

DO  2  J=2 , JMAX 

2  Z< J)=Z< J)/ZSCALE 
DYDX 1 ^—Z SCALE 
RETURN 

END 

SUBROUTINE  ENDS 

C  #*******#**#  »  »*************##*#*#***#####*****#'»'»*#•»***'**«  *«'#*******  • 
C  THIS  SUBROUTINE  COMPUTES  VELOCITIES  ON  THE  KEEL  AND  FREE  SURFACE 
C  *****************************  *******#**•*#*********•»(••*■*•*•*■*  *•*•*■*•**•***•»•■*•#■** 
COMMON/COORDB/XB ( 40 ) »  YB  <  40  >  i  ZB  <  40  ) 

COMMON/VEL/  UE < 40 > , VE < 40 ) » WE < 40 ) 

COMMON/NDATA/  IUli IU2» IU3» IU4» IU5, IMAX t JMAX . NXTLt NZTL t JFZ » LZ < 30 > 

Y2SQ=YB<2)**2 

Y3SO=YB<3>**2 

D-l. /(Y3S0-Y2SQ) 

UE  < 1 )  «  ( Y3SG*UE ( 2 ) -Y2SQ*UE ( 3 ) ) *D 
VE ( 1 >«0. 
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12100 
12200 
12300 
12400 
12500 
12600 
12700 
12800 
12900 
13000 
13100 
13200 
13300 
13400 
13500 
13600 
13700  C 
13800  C 
13900  C 
14000 
14100 
14200 
14300 
14400 
14500 
14600 
14700 
14800 

l^fVUVJ 

15000 
15100 
15200 
15300 
15400 
15500 
15600 
15700 
15800 
15900  2 

16000 
16100 
16200 
16300 
16400 
16500 
16600 
16700 
16800 
16900  3 

17000 
17100 
17200 
17300 
17400 
17500 
17600 
17700 
17800 
17900 
18000 


WE  <  1  )  =  <  Y3S0*WE  ( 2  )  - Y2SGHM4E  <  3  )  )  *D 
YB< 1 )~0. 

2B<1 )«<Y3S0*ZB(2)-Y2SQ*ZB<3>  )*D 
JM1  *=JMAX— 1 
JM2=JMAX— 2 
ZNM1  =  ZB ( JM1 ) 

ZNM2*=ZB  ( JM2 ) 

D»l. /(ZNM1-ZNM2) 

UE ( JMAX )  =  ( ZNM 1 *UE ( JM2 ) — ZNM2*UE  <  JM1 )  >*D 
VE ( JMAX )  =  ( ZNM1 *VE ( JM2 ) — ZNM2*VE  <  JM1 )  )*D 
WE ( JMAX ) =0, 

YB< JMAX )«< ZNM1*YB< JM2 ) -ZNM2*YB < JM1 ) >*D 
ZB< JMAX )=0. 

RETURN 

END 

SUBROUTINE  INITL ( X i Y» N, Li JF, ISTART, DYDX1 ) 
************************************************ ******************** 
THIS  SUBROUTINE  SETS  UP  INITIAL  VALUES  FOR  THE  VARIOUS  ARRAYS 
******************************************************************** 
DIMENSION  X (40) , Y<40) , L<30> 

COMMON /QU AN/  A  <  30  )  »  B  <  30 )  i  PH0PH0  (  30  )  *  PH0PH1  (  30  )  , 

*  PH0PS0  <  30 ) i PH0PS1 (30) »  PHI  PHI ( 30 >  » PH1PS0<30> , 

*  PH 1  PS 1 ( 30 ) i P 1 <  30 ) *  P2  <  30 ) ,  P3 ( 30 >  »  Q 1  < 30 ) , 

*  02(30) , 03(30) , XI (30) ,  YPHI0<30) , YPHI 1 (30) , 

*  H ( 30 ) 

COMMON/DIMN/  JFM, NFS 
COMMON/LDATA/  ISYMi ISET 
JFM-JF— 1 

Y1=Y<  X  ) 

YN-Y ( N ) 

IF ( I START- NE- 0 )  GO  TO  2 
X32=X  <  3 ) — X ( 2 ) 

X31=SX  <  3  )  — X  <  1  ) 

X21=X(2)-X( 1 ) 

AA=— ( X21  +  X31 ) / <  X21*X31 ) 

BB=X31/(X21*X32) 

CC*-X21/(X31*X32> 

DYDX 1 ~AA*Y ( 1 ) +BB*Y ( 2 ) +CC*Y ( 3 ) 

X32=X  <  N— 2 ) —X  <  N— 1 ) 

X31 =X  <  N— 2 )  —  X ( N ) 

X21—X (N— 1 )-X(N) 

AA~-(X21+X31 )/(X21*X31 ) 

BB*X31/ (  X21*X32> 

CC=— X21 / ( X3 1  *  X32 ) 

DYDXN=AA* Y ( N ) +BB*Y ( N- 1 ) +CC* Y ( N-2 ) 

IF ( I SET. EO. 0 )  GO  TO  3 
DYDXN=0. 

YN*— <  SB* Y ( N- 1  ) +  CC*Y  <  N-2  > ) /AA 
ISET*0 
PI < 1 >=Y1 
P2< 1 )«0. 

P3( 1 )*0. 

01 ( 1 >~DYDX1 
02  < 1 )«0. 

03 < 1 )«0. 

A  (  JF  )  **YN 
B  <  JF  )  *=DYDXN 
DO  1  J«liJF 
K»L< J) 

XI < J)«X<K) 


18100 
18200 
18300 
18400 
18500 
18600 
18700 
18800 
18900 
19000  1 

19100 
19200 
19300 
19400  C 
19500  C 
19600  C 
19700  C 
19800 
19900 
20000 
20100 
20200 
20300 
20400 
20500 
20600 
20700 
20800 
209uu  X 
21000 
21100 

21200 
21300 
21400 
21500 
21600 
21700 
21800 
21900 
22000 
22100 
22200  3 

22300 
22400 
22500 
22600 
22700  2 

22800 
22900 
23000 
23100 
23200 
23300 
23400 
23500 
23600 
23700 
23800 
23900  C 
24000  C 


PH0PH0 ( J ) =0 • 

PH0PH1 ( J>=0. 

PH0PS0 ( J ) =0 . 

PH0PS1 ( J>=0. 

PHI  PHI ( J>=0. 

PHI PS0  <  J ) =0. 

PH1PS1 ( J)=0. 

YPHI0  (J>=0. 

YPHI1  <J)=0. 

CONTINUE 

RETURN 

END 

SUBROUTINE  INTERP (IvJtXXt I FLAG* II) 

********************************************************************* 
THIS  SUBROUTINE  USES  THE  SPLINE  FIT  TO  INTERPOLATE  THE  ORIGINAL 
DATA  AND  TO  COMPUTE  ITS  FIRST  DERIVATIVE 
********************************************************************* 
COMMON /OUAN/  A(30)iB(30), PH0PH0<30>  * PH0PH1 ( 30 ) » 

*  PH0PS0 ( 30 ) * PH0PS1 (30)i PHI  PHI (30) , PH1PS0<30) i 

*  PH1PS1 (30) *  PI (30) i P2(30>  »  P3(30) *  01 (30) i 

*  02 ( 30 ) 5  03 ( 30 ) *  XI (30) » YPHI0(30) » YPHI1 ( 30 ) » 

*  H ( 30 ) 

COMMON/ DER I V/  DYBDZ (40) , DZBDZ (40) , DWDZ1 (40) »  DWDZN(40> 

COMMON/ VEL/  UE ( 40 ) *  VE ( 40 ) »  WE ( 40 ) 

COMMON/ VEL I /  UEI (40i 40 ) i WEI (40*  40) 

COMMON/D I MN/  JFM*  NFS 
IFCXX.LE. XX (1+1 ) )  GO  TO  1 
1  =  1  +  1 
X  M=X  + 1 

DX=(XX-XI ( I ) ) /H ( I ) 

DXP=( XI ( IP)-XX ) /H< I ) 

IF ( I FLAG. GT. 5 )  GO  TO  2 
PHI 1 = ( 3. -2. *DX ) *DX**2 
PHI0=1.-PHI1 
PSI 1=-H( I ) *DX  P*DX**2 
PS I 0=H ( I )*DX*DXP**2 

YY=A( I ) *PH I 0+ A ( IP)*PHI 1+B( I )*PSI0+B( IP)*PSI 1 
I F ( I FLAG • EO .  0 )  UEI ( 1 1  *  J ) =YY 
IF ( I FLAG. NE. 1 )  GO  TO  3 
WEI (I  I* J)=YY 

IF ( J. EO-  1 . OR. J. EO*  NFS  >  GO  TO  2 
IF ( I FLAG- EO- 2)  UE( J ) =YY 
I F ( I FLAG . EO • 3  >  WE(J)=YY 
IF ( I FLAG- EO. 4 )  DWDZ1 ( J ) =YY 
IF  < IFLAG. EO. 5 )  DWDZN( J)=YY 
RETURN 

PHI 1 P=6. *DX* ( 1  -  — DX ) /H ( I ) 

PHI0P=-PHI1P 
PSI1P=-DX*(2. *DXP— DX ) 

PSI0P=DXP*(DXP-2. *DX) 

YP=A< I ) *PHI0P+A (IP) *PHI 1 P+B  < I )*PSI0P+B< IP)*PSI1P 

IF < IFLAG. EO-  1 . AND. J-EO-  1 )  DWDZ1 ( I  I )=YP 

IF ( IFLAG. EO. 1 . AND. J. EO.NFS)  DWDZN( I I >=YP 

I F  < I FLAG . EO • 6  >  DYBDZ (J)=YP 

IF  < IFLAG. EO. 7 )  DZBDZ(J)«YP 

RETURN 

END 

SUBROUTINE  NAM IN 

******************************* ************************** ************ 
THIS  SUBROUTINE  COORDINATES  READING  OF  INPUT  DATA 
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24100  C  ********************************************************************  i 
24200  C0MM0N/NDATA/  N(40) 

24300  COMMON/MDATA/  M(31) 

24400  COMMON/LDATA/  ISYMt ISET 

24500  ISYM=0 

24600  CALL  NAMLST (6v2i I SYM , DUMMY ) 

24700  DO  1  1*1 '10  ~J 

24800  N<I)  =  0 

24900  1  CALL  NAMLST (6,2,N( I ), DUMMY ) 

25000  I  UP=  1 0+N (10) 

25100  DO  2  1  =  11, IUP 

25200  N  < I ) =0  ^ 

25300  2  CALL  NAMLST < 6, 2, N ( I ) , DUMMY ) 

25400  M(1>=0 

25500  CALL  NAMLST (6i2iM( 1 ) » DUMMY > 

25600  IUP= 1 +M ( 1 ) 

25700  DO  3  1=2, IUP 

25800  M ( I ) =0 

25900  3  CALL  NAMLST  (6,2,M(I),  DUMMY )  '{ 

26000  RETURN 

26100  END 

26200  SUBROUTINE  NAMLST (II,  10,  N,  X  ) 

26300  C  ******************************************************************** 
26400  C  THIS  SUBROUTINE  READS  INPUT  DATA 

26500  C  ******************************************************************** 
26600  DIMENSION  A (3) 

26700  IF ( N  .  NE.  0)  GO  TO  1 

26800  READ  <IIi2)  A(l)»A(2)iA(3)iN 

26900  WP  T  TF  <  10,2)  A(M»A(2)»A(3)iN 

27000  RETURN 

27100  1  READ  (11,3)  A< 1 ) t A ( 2) *  A (3 ) , X 

27200  WRITE  < 10,3)  A < 1 ) , A ( 2 > , A ( 3 ) *  X 

27300  RETURN 

27400  2  FORMAT ( IX, 3A4, 14) 

27500  3  FORMAT ( 1 X i 3A4  ,  1  PE 13.6) 

27600  END 

27700  SUBROUTINE  ROTATE (J) 

27800  C  ******************************************************************** 
27900  C  THIS  SUBROUTINE  TRANSFORMS  VELOCITY  COMPONENTS  TO  THE  NONORTHOGONAL 
28000  C  COORDINATE  SYSTEM 

28100  C  ******************************************************************** 
28200  COMMON/VEL/  UE ( 40 ) ,  VE ( 40 ) , WE ( 40 ) 

28300  COMMON/VEL I /  UEI (40,40) ,  WEI (40,40) 

28400  COMMON/NDATA/  IU1,  IU2,  IU3,  IU4,  IU5,  I MAX , JMAX  ,  NXTL , NZTL , JFZ  *  LZ  ( 30 ) 

28500  DO  2  1=1, NXTL 

28600  READ ( IU4)  HI , A, A , A, At  A, THETA, A 

28700  CSCTH= 1 • /S IN ( THETA ) 

28800  COTTH=CSCTH*COS ( THETA ) 

28900  UEI ( I, J)=UE( I ) —WE ( I )*COTTH 

29000  WEI ( I , J ) =WE ( I ) *CSCTH 

29100  IF ( J- EQ.  1. OR. J.EQ. NZTL)  WEI ( I , J  >  =0. 

29200  2  CONTINUE 

29300  RETURN 

29400  END 

29500  SUBROUTINE  SMOOTH ( X , Y, L , JF , N, ISTART , NN, XT, I FLAG, DYDX 1 , II) 

29600  C  ******************************************************************* » 
29700  C  THIS  SUBROUTINE  FITS  A  TABULATED  FUNCTION  WITH  CUBIC  SPLINES  USING 
29800  C  LEAST-SQUARES  ERROR  MINIMIZATION 

29900  C  ******************************************************************** 
30000  DIMENSION  X ( 40 ) , Y ( 40 ) , L ( 30 ) , XT ( 40 ) 
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30100  COMMON/QUAN/  A ( 30 ) . B ( 30 ) . PH0PH0 ( 30 ) * PH0PH1 ( 30 ) » 

30200  *  PH0PS0  <  30 ) »  PH0PS 1 ( 30 ) i PH 1  PH 1 <  30 ) i PH 1 PS0  <  30 ) » 

30300  *  PH1PS1  (30)  ,  PI  (30)  i  P2(30>  i  P3(30>  .  <51  (30)  » 

30400  *  <52(30)  .  (53(30)  i  XI  (30)  .  YPHI0(30)  »  YPHI 1  (30) » 

30500  *  H(30) 

30600  COMMON/D I MN/  JFM. NFS 

30700  C  INITIATE  THE  COMPUTATION 

30800  CALL  INITL (X.Y.N.L.  JF . I START  *  DYDX 1 ) 

30900  CALL  SUMUP(X» Y.N.L. JF) 

31000  C  COMPUTE  THE  P’S  AND  <5’S 
31100  DO  1  J=2»  JFM 

31200  M*J-1 

31300  C  COEFFICIENTS  IN  THE  ORIGINAL  EQUATIONS 
31400  A1=PH0PH1(M> 

31500  B1 =PH1 PHI ( M  >  +PH0PH0 ( J ) 

31600  C1=PH0PH1(J> 

31700  D1-PH1PS0(M> 

31800  E1=PH1 PS1 ( M ) +PH0PS0 ( J ) 

31900  F1=PH0PS1(J> 

32000  G1»YPHI1  ( M ) +YPH I 0  (J) 

32100  A2=3. /H(M>**2 

32200  B2=3 . * ( 1 . /H ( J  > **2- 1 . /H ( M ) **2 ) 

32300  C2~-3. /H( J>**2 

32400  D2=1./H(M) 

32500  E2=2. *( 1 . /H( J)+l . /H(M) ) 

32600  F2=1./H(J> 

32700  G2=0. 

32800  C  MODIFIED  COEFFICIENTS 
32900  6 1  b=&  1  -*•  A  1  *P2  (  M  )  +0 1  <  M  ) 

33000  E1S=E1+A 1 *P3 ( M ) +D 1 *<53 ( M ) 

33100  G1S=G1-A1*P1  (M)-D1*<51  (M) 

33200  B2S=B2+A2*P2 ( M  >  +D2*<52 ( M ) 

33300  E2S=E2+A2*P3 ( M ) +D2*<53 ( M ) 

33400  G2S=G2— A2*P 1  (M)-D2*<51  (M) 

33500  C  THE  P’S  AND  <5’S 

33600  DEL=1. / ( B 1 S*E2S— B2S*E IS) 

33700  PI ( J)«(E2S*G1S-E1S*G2S)*DEL 

33800  P2( J)*(E1S*C2  -E2S*C1  >  *DEL 

33900  P3( J)*(E1S*F2  -E2S*F1  )*DEL 

34000  <51  (  J)=(B1S*G2S-B2S*G1S)*DEL 

34100  <52(  J)»(B2S*C1  -B1S*C2  )*DEL 

34200  <53  ( J)«(B2S*F1  -B1S*F2  >*DEL 

34300  1  CONTINUE 

34400  C  SOLVE  FOR  THE  A’S  AND  B’S 
34500  K»JFM 

34600  KP*JF 

34700  DO  2  J— 1 . JFM 

34800  A ( K ) “PI (K)+P2(K)*A(KP)+P3(K)*e<KP) 

34900  BOO  =<51  00+<52<K>*A<KP>-*-<53(K>*B<KP> 

35000  KP=K 

35100  K=K-1 

35200  2  CONTINUE 

35300  C  INTERPOLATE  TO  DESIRED  MESH 
35400  1  =  1 

35500  NFS-NN 

35600  DO  3  J=1.NN 

35700  CALL  INTERP (IiJ»XT(J)» I FLAG i 1 1 ) 

35800  3  CONTINUE 

35900  RETURN 

36000  END 
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36100 
36200  C 
36300  C 
36400  C 
36500  C 
36600 
36700 

36800 

36900 

37000 

37100 

37200 

37300 

37400 

37500 

37600 

37700 

37800 

37900 

38000 

38100 

38200 

38300 

38400 

38500 

38600 

38700 

38800 

38*00 

39000 

39100 

39200 

39300 

39400 

39500 

39600 

39700 


SUBROUTINE  SUNUP < X , Y, Ni Li JF ) 

****************-**•****■*****■***•****■******•************■**  **-*-**■-**■***■*■-*•** 
THIS  SUBROUTINE  PERFORMS  SUMMING  OPERATIONS  ON  THE  CUBIC-SPLINE 
WEIGHTING  FUNCTIONS 


*#**********************-********************•***■**#*.*.***  ************* 
DIMENSION  X ( 40 )  »  Y  ( 40 ) i L ( 30 ) 

COMMON/OUAN/  A<30) , B<30> i PH0PH0<30>  *  PH0PH1 (30) i 


1 

2 


* 

* 


PH0PS0 ( 30 ) i PH0P SI (30) *  PH 1  PH 1 ( 30 ) i PH1PS0(30) t 
PH1PS1 <  30 ) i P 1 <  30 ) ?  P2 ( 30 ) »  P3  <  30 ) *  0 1 ( 30  > i 


■* 


Q2(30) ,03(30) i XI (30) , YPHI0C30) ,  YPHI1 (30) , 


* 


H<  30  ) 


COMMON/DIMN/  JFMi NFS 
DO  2  J= 1 i JFM 
H ( J ) =X I ( J+l )-XI ( J) 


KL=L< J) 

KU=L( J+l ) 

DO  1  K»KLiKU 


DX=( X (K)-XI ( J) ) /H( J) 
DXP=  (XI  ( J4- 1  )  —  X  (  K  )  )/H(J) 
YY=Y(K) 


PHI 1 = ( 3. -2. *DX )*DX**2 


PH 1 0= 1 . —PHI  1 

PSI 1®- H( J)*DXP*DX**2 


PSI0=  H(J 
PH0PH0 ( J ) 
PH0PH1 ( J) 
PH0PS0 ( J ) 
PH0PS1 ( J) 

PM  1  PHI  <  ..T  ) 

PHI PS0 ( J ) 
PHI PSI ( J) 
YPHI0  (J) 
Y  PH I  1  (J) 

CONTINUE 
CONTINUE 
RETURN 
END 


)*DX*DXP** 
=PH0PH0 ( J ) 
=PH0PH 1 ( J) 
=PH0PS0 ( J ) 
=PH0PS 1 ( J) 
=  PH1  PHI 
—  PH 1 PS0  <  J ) 
—PHI PSI ( J) 
=YPHI0  (J) 
=YPHI1  <J) 


+PH I 0*PH I 0 
+PHI0*PHI 1 
4  PH 1 0*-PS  1 0 
PH  1 0*  PS  I  1 

-4-  PM  T  1  ♦PMT  i 

4-PHI 1*PSI0 
+PHI1*PSI 1 
+  PH I 0* YY 
4-PHI  1*YY 


00100  C  PROGRAM  EDDY3 

00200  C  *********************************************  ***■#  *•*#•*«•*  **#*■*  -»*# ******* 
00300  C  THREE-DIMENSIONAL  BOUNDARY  LAYER  PROGRAM  FOR  FLOW  ABOUT  SHIP  HULLS 
00400  C  USING  THE  W I LCOX-RUBES I N  TWO-EQUATION  MODEL  OF  TURBULENCE 

00500  C  *********************•********************•*******■**********••*-***  v***** 
00600  INTEGER  TAPE I N i TAPEOT , TAPEGP , TAPEPF, TAPEDT , TAPEVL 

00700  C  DIMENSION  I NDEX 1(11 89 ) , I NDEX4 (11 89 ) 

00800  COMMON/ BLC0/  NX1NZ1NP1 NXT ,  NZT i NTRi NPT , NXTL , NZTL ,  NXSTRT , 

00900  *  NZSTRT,KC, IT, I FLOW, ICHORD, ISPAN, INTDIR 

01000  COMMON /BLC3/  X (01 ) , Z<41 ) , DETA ( 101 ) , ETA ( 101 ) , ETA E, VGP, CEL, BEL 1 , BEL 

01100  C0MM0N/BLC5 /  F(101,2,2),U<101,2,2),V(101),G<101,2,2),W(101,2,2), 

01200  *  T( 101 ) , TPROF ( 620 ) ,TPCF< 10) 

01300  COMMON /BLC9/  NSEP ( 4 1  ) ,  I  CASE ( 4 1 ) 

01400  COMMON/TURB/  E< 101 , 2,2) , WT ( 101 ) , WT2 ( 101 ,2,2), ELT ( 101 ) 

01500  COMMON/STOR/  CON, EPS, I CONVE » I TMAX , WTEDG 

01600  COMMON/TAPE/  TAPEIN, TAPEOT, TAPEGP, TAPEPF , TAPEDT , TAPEVL 

01700  COMMON/ I PRT /  IPZ,  IPX,  I  PR I NT , EPSV , EPST 

01800  COMMON/ CONV/  ICOUNE, ICOUNWt INE6E, INEGW, DELV, DELT 

01900  COMMON /SVNM/  IPCF,IPRF 

02000  DATA  IPCF, IPRF/10, 620/ 

02100  C  CALL  OPENMS< 16, INDEX1 , 1 189, 0) 

02200  C  CALL  OPENMS( 17, INDEX4, 1 189, 0 ) 

02300  C*-**  DISK  FILES  TAPEGP(16)  AND  TAPEPF<17)  ARE  RANDOM  ACCESS  FILES 
02400  C  TAPEGP. .. BODY  GEOMETRY  AND  VELOCITY  DISTRIBUTION 

02500  C  TAPEPF. .  . PROFILES  AND  PRESSURES 

02600  C***  DISK  FILE  TAPEDT (18)  IS  CREATED  IN  PROGRAM  SHPMSH 

02700  C***  DISK  FILE  TAPEVLC19)  IS  CREATED  IN  PROGRAM  VELOC 

02800  C***  FILES  TAPE  I N ( 5 )  AND  TAPE0T<6>  ARE  CONVENTIONAL  READ  AND  WRITE 

0^vtow  u  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  ~ 

03000  CALL  INPUT 

03100  CALL  GRID 

03200  DO  10  I - 1 , NZTL 

03300  10  NSEP ( I ) =0 

03400  DO  20  K«l,2 

03503  DO  20  J= 1 , NPT 

03600  F( J«K« 1 )*0. 

03700  U( J,K, 1 >=0. 

03800  G< J, K, 1 )^0. 

03900  “  W(J,K,1)=0. 

04000  E <  J , K,  1  ) -0 - 

04100  20  WT2< J,K, 1 )=0. 

04200  I NTD I R~  I FLOW 

04300  IF  < I FLOW. EQ.  1  )  WR I TE ( TAPEOT , 690 ) 

04400  IF  < I FLOW • EQ. 2)  WR I TE ( TAPEOT , 680) 

04500  NZ 1 =NZSTRT 

04600  NZ2=NZT 

04700  DO  500  NX=NXSTRT , NXT 

04800  NZ-NZ 1 

04900  I F  < I FLOW • EQ . 2 )  NZ=NZ2 

05000  CALL  LOG  I C  <  NZ 1 , NZ2 ) 

05100  30  IF(NX,EQ. NXSTRT)  CALL  PROFIL(l) 

05200  IF  <  NSEP  <  NZ ) . NE . 0 )  GO  TO  405 

05300  I T =0 

05400  I CONVE- 0 

05500  I  PR I NT  =  0 

05600  IF< (NZ. EQ. NZSTRT. OR. NZ. EQ. NZT ) . AND. MOD (NX, IPX ) . EQ. 0)  IPRINT=1 

05700  I F  <  NX . EQ . NX T . AND . MOD ( NZ ,  I P Z ) . EQ . 0 )  I  PR  I NT  =  1 

05800  IF  <  MOD  (  NZ  ,  I  PZ  )  .  EQ.  0.  AND.  MOD  ( NX  ,  IPX).EQ.0)  IPRINT*=1 

03900  KC-2 

06000  IF  < I CASE ( NZ ) . EQ. 4 )  KC  =  4 
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IF  <  NZ .  E<5.  NZT  )  KC=ICHORD 

IF<NZ.  EQ.  NZ8TRT)  KC=1 

IF < I PRINT. EQ. 0 )  GO  TO  290 

WRITE  <  T  APEOT i 620 )  NXiNZi X  ( NX  ) i Z(NZ) 

IFCKC.EQ.  1  )  WRITE  <  T  APEOT  >  640) 

IF  <  KC- E0.2)  WR I TE ( T APEOT i 650 ) 

IF < KC. EQ. 4 )  WR I TE  <  T APEOT i 670 ) 

290  IT=IT+1 

IF< IT.LE. ITMAX)  GO  TO  300 
WRITE ( T APEOT  %  630 )  NX*NZ 
GO  TO  390 
300  CALL  EDDY 
CALL  CGEN 
CALL  S0LV6 
CALL  SOLVEW 

IF  <  V  < 1 ) . LT -  0-  )  GO  TO  380 

IF <  ABS ( DELV/ <  V ( 1 )+.5*DELU> ) . GT.EPSV.OR. 

*  ABS(DELT) . GT. EPST)  GO  TO  290 
IF < ICONVE. EQ. 0)  GO  TO  290 

IF ( NP. EQ. NPT )  GO  TO  390 

IF  ( U  <  NP— 1*2*2). LT -  . 999 )  GO  TO  370 

GO  TO  390 

370  CALL  PROF I L ( 3 ) 

I T =0 

GO  VO  290 

380  WRITE ( TAPEOT i 720 )  NXiNZ 
CALL  SHIFT2 <  NZ 1 »  NZ2 ) 

GO  TO  405 

37S  CALL  OUTPUT ( i ) 

IFCNP- EQ. NPT)  GO  TO  400 
CALL  PROF I L  <  2 ) 

400  CALL  SHIFT <  NZ  1 i NZ2  > 

4F5  I F  < I FLOW • EQ •  2 )  GO  TO  407 
NZ=NZ+1 

IF<NZ.GT.NZ2>  GO  TO  410 
GO  TO  30 
407  NZ=NZ— 1 

IF ( NZ • LT • NZ 1 )  GO  TO  410 
GO  TO  30 

410  CALL  OUTPUT  <  2 ) 

500  CONTINUE 

620  FORMAT  <  1H0*  5H*--**-*  i  4HNX  =iI3*5X4HNZ  =iI3i5X3HX  F10.  3»  5  X3HZ  =  * 

*  F 1 0 , 3  ) 

630  FORMAT ( 1H0i 1 8X32H I TERAT I ONS  EXCEED  ITMAX  FOR  NX  =iI4i5X4HNZ  =i I 
640  FORMAT < 1H0i 30H**  SYMMETRY-PLANE  EQUATIONS  **> 

650  FORMAT  <  1H0*45H**  CHORDWISE  I NF I N I TE-SWEPT  HULL  EQUATIONS  ) 

670  FORMAT  (  1H0*  1QH**  GENERAL  CASE  -**  ) 

680  FORMAT( //1H0,  38X50H*-**** - WATERLINE  CALCULATIONS  STARTED - 

*■***//  ) 

690  FORMAT ( //1H0> 38X52H***** -  KEEL-LINE  CALCULATIONS  STARTED  — 

**•***■*//  ) 

720  FORMAT ( 1H0* 1 1X36H**  BOUNDARY-LAYER  SEPARATION  AT  NX  «»I4,5X4HNZ 
•*  T  4 ) 


00100  SUBROUTINE  ARRYMV < LL *  I , K ) 

00200  COMMON /BLC0/  NX , NZ , NP, NXT, NZT,  NTR,  NPTi  NXTl.»  NZTL*  NXSTRT , 

00300  *  NZSTRT.KC, IT, IFLOW, ICHORD, ISPAN, INTDIR 

00400  COMMON/BLC5/  F < 10 1 ,2,  2 > , U < 10 1 ,2,  2 ) , V ( 10 1 > , G < 10 1 , 2,  2 ) , W < 101 » 2, 2 )  , 

00500  *  T ( 1 01 ) , TPROF ( 620 ) . TPCF (10) 

00600  COMMON/PRSC/  PS 1 ( 2,  2 > , HI ( 2 , 2 > , H2 < 2. 2 ) , UE < 2, 2 > , 

00700  *  WE (2.2) . CK1 (2,2). CK2 (2.2) . CK 12(2,2) . CK21 (2.2). 

00800  *  THETA (2.2), PFRS , UFRS »  CNUFRS, UREF 1 , WNP 

00900  COMMON/PRES/  PI (2.2), P2 (2,2) , P3 (2.2) , P4<2, 2) , P5  <  2 ,  2  ) , 

01000  *  P6<2»2> , P7(2, 2) » P8<2,2) , P9<2, 2) » P1012, 2) , PI  1 (2. 2) , 

01100  *  P1212, 2) , P13<2, 2) , P14<2, 2) 

01200  COMMON/TURB/  E < 101 , 2, 2 ) , WT < 10 1 ) , WT2 < 10 1 , 2, 2 > , ELT < 101 ) 

01300  COMMON/UN I V/  J1 ,  J2,  J3,  J4,  J5  ,  N1 ,  N2»  N3,  N4,  N5  ,  N6,  N7,  N8,  N9, 

01400  *  N10, N1 1 , N12, N13, N14, NPT2, NPT3, NPT4, NPT5, NPT6 

01500  C---------------------------------- 

01600  NPT2=NPT*2 

01700  NPT3=NPT*3 

01800  NPT  4=NPT*4 

01900  NPT5- NPT  *5 

02000  NPT  6=NPT  *6 

02100  GO  TO( 10, 30, 90, 50, 70, 85 > ,  LL 

02200  C  PUT  FLOW  ARRAYS  IN  TPROF  IN  PREPARATION  FOR  A  DISK  WRITE 
02300  10  DO  20  J~ 1 , NPT 

02400  CALL  INDEX (J, NPT) 

02500  TPROF ( J  >  =F ( J ,  I , K ) 

02600  TPROF ( J 1 ) =U ( J »  I , K  > 

02700  TPROF ( J2 ) — G  <  J,  I , K ) 

02800  TPROF < J3)=W< J, I.K) 

02900  T PROP  tJ'*)  — 1. 1  J  ,  1 ,  K  > 

03000  20  TPROF ( J5 ) =WT2 ( J, I.K) 

03100  TPROF INI )=P1 ( I.K) 

03200  TPROF  <  N2  )  =  P2  (I.K) 

03300  TPROF CN3>~P3< I.K) 

03400  TPROF <N4)=P4( I.K) 

03500  TPROF ( N5 ) — P5 ( I , K ) 

03600  TPROF ( N6 ) =P6 (I.K) 

03700  TPROF ( N7 ) =P7 ( I , K ) 

03800  TPROF (N8)=P8( I.K) 

03900  TPROF <N9)~P9f I , K) 

04000  TPROF(N10)=P10< I ,K) 

04100  TPROF  (Nil  )==P11<  I,K) 

04200  T  PROF ( N 1 2 ) - P 1 2 ( I.K) 

04300  TPROF (N13)=P13( I , K ) 

04400  TPROF (N14)~P14( I.K) 

04500  RETURN 

04600  C  FILL  FLOW  ARRAYS  AFTER  A  DISK  READ 
04700  30  DO  40  J=1 , NPT 

04800  CALL  INDEX (J, NPT) 

04900  F(J, I,K) = TPROF ( J ) 

05000  U< Ji I,K)*TPROF< J1 ) 

05100  G(J, I,K)»TPPOF( J2) 

05200  W( J, I ,K)=TPROF( J3) 

05300  E ( J, I.K) =TPROF ( J4 ) 

05400  40  WT2 ( J , I.K) =TPROF ( J5 ) 

05500  PI  (  I ,  K)«=TPROF<Nl  ) 

05600  P2 (I.K) =TPROF ( N2 ) 

05700  P3 (I.K) =TPROF  <  N3 ) 

05800  P4 ( I , K ) =TPROF ( N4 ) 

05900  P5 ( I , K ) *TPROF ( N5 ) 

06000  P6< I,K)=TPR0F(N6) 
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P7 ( X  v  K  > “TPROF ( N7 ) 

P8< 1 1 K  > =TPROF  <  N8 ) 

P9 (IiK) =TPROF <N9) 

P10< IiK)=TPROF(N10) 

Pll(IiK) =TPROF (Nil) 

P12< IiK) =TPROF ( N 1 2 ) 

P13(IiK) =TPROF ( N 1 3 ) 

PI 4 ( I i K ) =TPROF ( N 1 4 ) 

RETURN 

C  SHIFT  (K?K)  FLOW  ARRAYS  TO  (IiK) 

50  DO  60  J=1 i NPT 

F( Ji  IiK) =F ( J1K1K) 

U< Ji I i K)«U( Ji Ki K) 

G (Ji IiK)=G( JiKiK) 

W ( Ji IiK)»W(JiK»K) 

E  ( Ji  1 1 K ) =E (JtKiK) 

60  WT2( Ji I i K)-WT2< Ji Ki K) 

RETURN 

C  SHIFT  ( I i K )  FLOW  ARRAYS  TO  <K»K> 

70  DO  80  J= 1 i NPT 

F< JtKiK) =F <  J i  1 1 K ) 

U( Ji Ki K)=U( Ji  1 1  K ) 

6(JiKiK)=6(Ji  1 1 K  > 

W< JiKiK)«W( Ji IiK) 

E( J'KiK)«E< Jf IiK) 

80  WT2< JiKiK)«WT2( Ji IiK) 

RETURN 

C  SHIFT  <IiI)  FLOW  ARRAYS  TO  (IiK) 

85  00  88  J«1,NPT 

F< Ji I i K)=F< Ji Iil) 

U< Ji IiK) =U ( J i Iil) 

G< Ji IiK) =6 ( J i Iil) 

W( Ji IiK) =W ( J i Iil) 

E< Ji IiK)=E<Ji  Iil) 

88  WT2( Ji I i K)=WT2< Ji I i I ) 

RETURN 

C  FILL  GEOMETRY  ARRAYS  AFTER  A  DISK  READ 
90  HI ( IiK)=TPCF< 1 ) 

H2< IiK) =TPCF ( 2 ) 

CK1 ( IiK) =TPCF ( 3 ) 

CK2 ( IiK)=TPCF(4> 

CK12 (IiK) =TPCF ( 5 ) 

CK21 ( I i K ) —TPCF ( 6 ) 

THETA ( IiK) =TPCF ( 7 ) 

PS1 ( IiK)=TPCF(Q) 

UE( IiK) =TPCF ( 9 ) 

WE ( I i K)=TPCF( 10) 

RETURN 

END 


4 


1 
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10 

20 


22 


25 
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30 


32 


35 

37 

40 


SUBROUTINE  BELCEL 

COMMON / BL C0 /  NX , NZ * NP , NXT , NZT , NTR , NPT , NXTL , NZTL , NXSTRT , 
*  NZSTRT t  KC ,  IT,  IFLOW,  I  CHORD,  ISPAN,  INTDIR 


C0MM0N/BLC3/ 
COMMON/BLC5 / 

* 

COMMON/PRES/ 

* 

* 


X ( 81 ) , Z  <41 ) , DETA  < 101 ) i ETA < 101  >  , ETAE, VGP, CEL, BEL 1  *  BEL 
F< 101 , 2,2) ,U< 101 , 2,2) , V< 101 > , G< 101 ,2, 2) , W( 101 , 2,2) , 
T<  101  >  , TPROF  <  620 ) , TPCF  < 10) 

PI (2«2) «  P2(2f 2) i P3<2»2>  f  P4<2v  2) v  P5(2«  2) * 

P6<2»  2) i P7(2t 2) i P8(2f  2) »  P9<2t  2) f  P10(2i 2> »  PI  1 <2i 2) i 
P12(2,2>,P13<2,2),P14<2,2> 


CALL  PRESUR 

IF< KC. EQ. 4 >  GO  TO  20 

BEL 1=0. 

BEL2=0. 

IF ( NX . EO- NXSTRT )  GO  TO  10 
CEL=P10(2,2)/(X(NX)-X<NX-1 ) ) 

RETURN 

CEL—P 1 0(2,2) /(X(NX+1 )~X(NX> ) 

IF ( NTR. EO. 0 )  CEL=0 . 

RETURN 

DEL  X  ~  X (NX)— X (  N  X  —  1 ) 

IF < I FLOW, EO. 2 )  GO  TO  30 
DEL  Z1=Z(NZ)~Z(NZ— 1 ) 

DELZ2—Z (NZ  +  1 )  ~Z  ( NZ ) 

IF  <  NZ .  GT  •  (  NZSTRT+  1  )  )  GO  TO  25 
DO  22  J  =  1 , NPT 
G  <  J,  2,  1  >*=0. 

W< J, 2, 1 )=0. 

GO  TO  40 

IF< ICH0RD.E0.2.0R.NZ.LT.  (NZT-1) >  GO  TO  40 
DO  27  J»1 «  NPT 
G< J, 1 i 1 >=0. 

W< Jt 1 i 1 )*0. 

GO  TO  40 

DELZ 1  =  Z (NZ )  —  Z ( NZ  + 1 ) 

DELZ2=Z ( NZ— 1 >-Z (NZ ) 

IF<  I  CHORD.  Et>i.  2.  OR.  NZ.  LT •  (NZT—l  )  )  GO  TO  35 
DO  32  J=  1 , NPT 
G< J, 2, 1 >=0. 

W< J, 2, 1 >=0. 

GO  TO  40 

IF(NZ.GT.  (NZSTRT  + 1  )  )  GO  TO  40 
DO  37  J= 1 ,  NPT 
G  ( J, 1 , 1 )=0. 

W<  J, 1 , 1 )=0. 

CEL«P10<2, 2) /DELX 
BEL1=. 5*P7<2, 2) /DELZ1 
BEL2- - 5*P7 (2,2) /DELZ 2 
RETURN 
END 
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SUBROUTINE  CGEN 


COMMON/ BLC0/ 
v 

C0MM0N/BLC3/ 
C0MM0N/BLC5 / 

* 

COMMON /PRSC/ 


* 

* 


50 


C0MM0N/BLC7/ 
COMMON/ PRES/ 


* 

* 


COMMON/BLCE/ 

COMMON/BLCY/ 


NX*NZ*  NP* NXT * NZT * NTR* NPT *  NXTL *  NZTL i NXSTRT * 

NZSTRT *  KCi IT* IFLOW* I  CHORD*  ISPAN*  INTDIR 
X ( 81 >  *  Z ( 4 1 ) * DETA ( 101 ) *  ETA ( 101 ) *  ETAE* VGP*  CEL*  BEL 1  *  BEL 
F( 101  *  2*  2) *U< 101  *  2*  2) *  V< 101 ) *  G< 10  *  2*2) * W< 101  *  2*2) * 
T< 101 >  * TPROF <  620 ) *  TPCF( 10) 

PS1 (2*2) *H1 (2*2) *  H2<2*  2) *UE(2*  2) * 

WE <2*  2) *  CK1 <2*  2) * CK2<2*2> *  CK12C2*  2) * CK21 (2*  2) » 

THETA  < 2*2) »  PFRS*  UFRS*  CNUFRS*  UREF1 *  WNP 
PU< 101 ) i PW< 101 ) *QU< 101 ) *GW< 101 ) 

PI (2*2) * P2<2* 2) *  P3<2* 2) *  P4(2*2) *  P5(2* 2) * 

P6(2*2) 9  P7(2*2) *  P8(2»2) *  P9(2* 2) 9  P10<2*  2) *  Pll (29  2) 9 
P12<2*2) 9  P13<2*2>  9  PI 4 <2* 2) 

EDV ( 101 ) 

Y1 ( 101 ) * Y3< 101 > * Wl* W2* W3* W4 


IF(IT.EQ.l)  CALL  BELCEL 
U ( NP*  2*  2  )  ~  1 « 

W ( NP *  2*2) =WE ( 2 *  2) /UREF1 
IF(KC.EO.l)  W ( NP*  2*  2 ) =WNP 
PU( 1 )=0. 

PW( 1 )=0. 

QU  (  1  )  =0 . 

OW( 1 >=0. 

NPM=NP- 1 

BM*1 . 5*<EDV( 1 >+EDV<2) ) 

DO  50  J=2 *  NPM 

BP= 1 . +. 5*  <  EDV ( J ) +EDV ( J+ 1 )  ) 
ijn  —  ij  (  J  -  2:2) 

WB~W ( J*  2*  2 ) 

E2=P 1 ( 2  *  2 ) *F ( J  *  2 *  2 ) + P6  <  2  *  2 ) *G  <  J  *  2  *  2 ) 

DFB=  CEL* ( F ( J  *  2  *  2 ) — F ( J  *  1*2)) 

DGe=BELl*(G< J* 2* 2)-G< J* 2*  1 )  ) +BEL2* ( G ( J  *  1  *  1 ) -G ( J *  1  *  2 )  ) 

USB=UB*UB 

WSB=WB*WB 

DVB= ( E2+DFB+DGB ) / ( DETA ( J ) +DETA ( J- 1 ) ) 

X 1 =- <  P2 ( 2 *  2 ) *UB+P5 (2* 2)*WB) 

X2=— ( P4 < 2* 2)*UB+P3<2*  2)*WB> 

A 1 =  BM*  Y3 ( J ) —DVB 

B 1 =-BP*Y 1 ( J ) -BM* Y3 ( J ) +X 1 - ( CEL*UB+BEL 1 *WB ) 

Cl=  BP*Y1 < J)+DVB 

D 1 =P8 ( 2  *  2 ) *WSB-P 1 1 <  2 *  2 ) - < CEL*UB+BEL2*WB ) *U  <  J  *  1  *  2  > 

*  —BEL  1  *WB*U  <  J  *  2  *  1  )  +BEL2*WB*U  (  J  *  1*1) 

B2-— BP*Y  1  <  J  )  -BM* Y3  <  J  )  +X2-  (  CEL*UB+BEL.  1  *WB  ) 

D2~P9 (2*2) *USB~  P 1 2 ( 2  *  2 ) - <  CEL*UB+BEL2*WB ) *W  <  J  *  1  *  2 ) 

*  —BEL 1 *WB*W ( J  *  2* 1 >+BEL2*WB*W< J* 1*1) 
B1S=B1+A1*QU< J-l ) 

B2S-B2+A 1 *OW <  J- 1 ) 

D 1 S=D 1 —A 1 *PU ( J— 1 ) 

D2S-D2-A 1 *PW  <  J- 1 ) 

PU( J)=D1S/B1S 

PW( J)=D2S/B2S 

QU< J)»-C1/B1S 

QW< J)»-C1/B2S 

BM-BP 

RETURN 

END 
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SUBROUTINE  EDDY 

COMMON /BLC0/  NX, NZ , NP, NXT, NZT, NTR, NPT, NXTL, NZTL, NXSTRT, 
*  NZSTRT i KC, IT, IFLOW, ICHORD, ISPAN, INTDIR 


COMMON/BLC5 / 


F< 101,2, 2), U( 101,2,2) ,  V  ( 101 ) ,G< 101, 2,2), W< 101,2,2), 


*  T< 101 > ,TPROF<620> ,TPCF< 10 > 

COMMON/ PRSC/  PS1 < 2, 2 ) , HI ( 2, 2 ) , H2 ( 2, 2 ) , UE ( 2, 2 > , 


* 


WE (2,2) , CK 1 <  2 , 2 ) , CK2 


(2,2) , CK1 2(2,2) , CK21 (2,2)  , 


COMMON /BLCE/ 
COMMON/TURB/ 
COMMON/STOR/ 


THETA <2, 2) , PFRS, UFRS, CNUFRS, UREF1 , WNP 
EDV( 101 ) 

E( 101,2,2) , WT( 101 ) , WT2 ( 101,2,2)  ,  ELT ( 101 ) 
CON, EPS, ICONVEt I TMAX , WTEDG 


REY-UE (2,2) *PS1 (2,2) /CNUFRS 
CEP— RE Y* ( UREF 1 /UE (2,2) ) 

EDU ( 1 )=0. 

ELT ( 1  >==0. 

DO  54  J =2 , NP 

IF ( IT. EQ. 1 )  WT ( J ) -SORT ( WT2 ( J , 2 , 2 ) ) 
RET=CEP*E <  J , 2 , 2 ) /WT ( J ) 

EDV ( J ) =  RET  * ( 1 . —CONVEX  P ( —RET ) ) 

ELT( J)=SORT(E( J, 2,2)) /WT ( J ) 

RETURN 

END 


00100 

00200 

00300 

00400 

00300 

00600 

00700 

00800  C  - 
00900 
01000 
01 100 
01200 
01300 
01400  30 

01500 
01600 
01700 
01800 
01900  40 

02000 
02100 
02200 
02300 
02400 
02500 
02600 
02700 
02800  50 

C  /-» 


SUBROUTINE  GRID 

INTEGER  TAPEINt  TAPEOT i TAPEGP»  TAPEPF*  TAPEDTi TAPEVL 
COMMON/BLC0/  NX  »  NZ  t NP»  NXT *  NZT i NTR, NPTi  NXTLi  NZTLi NXSTRTt 
*  NZSTRT »  KC» IT t IFLOW* ICHORDi ISPANi  INTDIR 

C0MM0N/BLC3/  X<81)iZ(41)« DETA< 101 >  «  ETA< 101 ) * ETAE»  VGPi CEL»  BEL1 ?  BEL 
COMMON/BLCY /  Y1 ( 101 > , Y3 ( 101 ) , W1 , W2, W3, W4  ^ 

COMMON/TAPE/  TAPEIN*  TAPEOT »  TAPEGPi TAPEPF*  TAPEDTt  TAPEVL 

ETA ( 1 ) *0.  ----- 

NP«=ALOG<  (ETAE/DETA(1)  )*(VGP-1.  >  +  i.  )  / ALOG ( VGP )  + 1 . 000 1 
IF(NP.LE.NPT)  GO  TO  30 
WRITE (TAPEOT »  50 > 

STOP  10 

DO  40  J=2»  NPT 

DETA ( J ) =VGP*DETA ( J- 1 ) 

ETA( J)=ETA( J-l )+DETA( J-l ) 

Cl =2. / ( DETA ( J ) +DETA  <  J— 1 ) ) 

Y1 ( J ) —Cl /DETA ( J ) 

Y3( J)=C1/DETA( J-l ) 

ETAE=ETA(NP> 

W5«  < 1 . +VGP ) * ( 1 . +VGP* ( 1 . +VGP ) ) *VGP**3 
DEN=1. /<W5*DETA< 1 > ) 

W4=< 1 . +VGP)*DEN 

W3= ( 1 . +VGP* ( 1 . +VGP ) ) **2*DEN 

W2=VGP* < 1 . + VGP ) *W3 

Wl  =  ( 1 . +VGP+ ( 1 . +VGP* ( 1 . +VGP) )**2*< VGP*< 1 . +VGP1-1 .  ) )*DEN 
RETURN 

(  1H0»  37HNP  EXCEEDED  NPT  —  PROGRAM  TERMINATED) 
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00100 

00200 

00300 

00400  C 

00500 

00600 

00700 

00800 

00900 

01000 

01100 

01200 

01300 

01400 

01500 

01600 

01700 

01800 

01900 

02000 

02100 

02200 

02300 

02400 

02500 

02600 


SUBROUTINE  INDEX ( J*  NPT ) 

COMMON/UN IV/  J 1 i J2* J3, J4» J5*N1 * N2» N3» N4, N5 t N6* N7, N8* N9, 

*  N10fNlliN12vN13iN14f  NPT2*  NPT3i NPT4»  NPT5 * NPT6 


Jl-J+NPT 

J2-J+NPT2 

J3«J+NPT3 

J4»J+NPT4 

J5-J+NPT5 

IF( J.LT.NPT)  RETURN 

N1-NPT6+1 

N2*N1+1 

N3—N2+1 

N4=N3+1 

N5—N4+ 1 

N6=N5+1 

N7~N6+1 

N8=N7+1 

N9-N8+1 

N10=N9+1 

Nil “N10+ 1 

N12-N1 1+1 

N13=N12+ 1 

N 1 4S=N1 3+  1 

RETURN 

END 


00100  SUBROUTINE  INPUT 

00200  DIMENSION  GMPKG ( 8 ) , TITLE ( 9 ) , APC ( 10 > 

00300  INTEGER  TAPE  IN. TAPEOT, TAPEGP, TAPEPF, TAPEDT, TAPEVL 

00400  COMMON/BLC0/  NX»NZ,NP*NXT,NZT»NTR,  NPT ,  NXTL, NZTL,  NXSTRT  » 

00500  *  NZSTRT,KC, IT, IFLOW, ICHORD, ISPAN, INTDIR 

00600  C0MM0N/BLC3/  X(81)»Z(41). DETA (101). ETA (101). ETAE ,  VGP ,  CEL , BEL 1 »  BEL 

00700  COMMON/ CHRD/  WNP1 ( 81 > , WNPN ( 81 ) 

00800  COMMON/ PRSC/  PS l(2»2)»Ht(2,2),  H2 (2,2) i UE (2,2), 

00900  *  WE(2, 2) , CK1 (2,2) , CK2(2, 2) , CK12(2,2) , CK21 (2, 2) » 

01000  *  THETA ( 2*  2  >  »  PFRS*  UFRS,  CNUFRS,  UREF 1 , WNP 

01100  COMMON/TAPE/  TAPEIN, TAPEOT, TAPEGP, TAPEPFi TAPEDT, TAPEVL 

01200  COMMON/ 1 PRT /  IPZ, IPX , I PRINT ,  EPS Vi EPST 

01300  COMMON/STOR/  CON* EPS, I CONVE, I TMAX , WTEDG 

01400  COMMON/ CONS /  S ALFA  *  SALF AS  *  SBET A »  SBET AS . SL AMD A »  SS I GMA  *  SSGMAS  , USTOP 

01500  *  ZIOTAE*  ZIOTAL.*  RW2 

01600  COMMON/SVNM/  IPCF*IPRF 

01700  COMMON/SHRT /  ISHORT 

01800  COMMON/ RELX/  RFTRB, RFVEL 

01900  COMMON/ INTG/  ALAMI ( 4 1 ) * CFX I ( 41 ) * CFZ I < 41 ) * DELTAI < 4 1 > * 

02000  *  THETXI (41 ) ,THETZI (41 ) 

02100  EQUIVALENCE  ( APC( 1 ) *  GMPKG ( 1 ) ) *  ( APC ( 9 ) * UEUF ) ,  ( APC ( 10 > , WEUF ) 

02200  DATA  TAPE IN /5 / *  TAPEOT /6/ *  TAPEGP/ 1 6/ *  TAPEPF / 1 7/ »  TAPEDT / 18/ * 

02300  *  TAPEVL/ 1 9/ *  I SHORT /0/ *  RFTRB/ . 8/ *  RFVEL/ 1 . / *  NPT / 101 / 

02400  DATA  SALFA/1 . 11111/*  SALFAS/. 3/*  SBET A/.  15/*  SBET AS/. 09/* 

02500  *  SLAMDA/ . 091 / »  SSI GMA/ . 5/ *  SSGMAS/ . 5/ *  USTOP/ 4. / * 

02600  *  ZIOTAE/3. 75E-5/, ZIOTAL/. 09/, RW2/4. / 

02700  DATA  EPS, EPSV, EPST /3* . 01 / *  ETAE/ 8. / *  I  CHORD/ 1 / , I FLOW/ 1 / , 

02800  *  IFPRNT/0/, ISPAN/1/* I TMAX /20/ * NTR/ 1 / 


02930 

NAriELIS  i  /NAht/ 

u£  l  A 

* 

fc.KB 

t 

EPST 

* 

ERBV 

t 

ETAE 

« 

03000 

* 

I  CHORD 

9 

IFLOW 

9 

IFPRNT 

f 

IPX 

f 

IPZ 

9 

03100 

* 

ISPAN 

9 

ITMAX 

9 

NPT 

9 

NTR 

9 

NXSTRT 

9 

03200 

* 

NXT 

9 

NZSTRT 

9 

NZT 

» 

TAPEIN 

9 

TAPEOT 

9 

03300 

* 

TAPEGP 

* 

TAPEPF 

i 

TAPEDT 

» 

TAPEVL 

9 

VGP 

9 

03400 

* 

ISHORT 

i 

RFTRB 

9 

RFVEL 

03500 

NAMELIST/DATA/ 

SALFA 

t 

SALFAS 

V 

SBETA 

SBETAS 

9 

SLAMDA 

9 

03600 

» 

SSI GMA 

9 

SSGMAS 

1 

USTOP 

* 

ZIOTAE 

9 

ZIOTAL 

1 

03700 

* 

RW2 

03800 

NAMEL I ST /STRT / 

ALAMI 

9 

CFX  I 

* 

CFZ  I 

9 

DELTAI 

9 

THETXI 

9 

03900 

* 

THETZI 

04000  c-------  -  __________ 

04100  READ(TAPEIN* 120)  TITLE 

04200  READ (TAPE  IN, NAME) 

04300  READ( TAPE IN, DATA) 

04400  READ ( TAPEIN, STRT ) 

04500  READ(TAPEIN, 170)  UFRS, PFRS* CNUFRS 

04600  READ (TAPEDT)  NXTL , NZTL 

04700  READ( TAPEDT)  ( X ( I ) , 1=1 , NXTL ) 

04800  READ ( TAPEDT )  ( Z ( I ) *  1  =  1  * NZTL ) 

04900  CON= 1 . -SLAMDA**2 

05000  IF( IFPRNT. EQ. 1 )  WRITE ( TAPEOT, 210 ) 

05100  NX=NXSTRT 

05200  IF( IFLOW. EQ. 1 )  NZ=NZSTRT 

05300  IF( IFLOW. EQ. 2)  NZ=NZT 

05400  UREF 1 =UFRS 

05500  I PNTG=8 

05600  DO  70  K= 1 , NZTL 

05700  DO  60  1=1, NXTL 

05800  IPNTG=IPNTG+1 

05900  READ (TAPEDT)  GMPKG 

06000  READ (TAPEVL)  UEUF, WEUF 


72 


06100  PE=PFRS 

06200  UEUF=UEUF*UFRS 

06300  WEUF=WEUF*UFRS 

06400  IF( IFPRNT. EQ. 1 )  WR I TE < TAPEOT . 220 )  Ki It GMPKG 

06500  CALL  STORIT ( TAPEGP*  APC* IPCF , I PNTGi  0 ) 

06600  IF ( I . NE. NXSTRT )  GO  TO  60 

06700  IF(K.NE.NZ)  GO  TO  60 

06800  HI (2,2>=APC( 1 ) 

06900  H2 ( 2>  2 ) =APC  <  2 ) 

07000  CK1 (2,2>=APC(3> 

07100  CK2 ( 2t  2) =APC ( 4 ) 

07200  CK12  (  2,2)  =APC  ( 5  ) 

07300  CK21 (2,2>=APC(6> 

07400  THETA ( 2*  2 ) =APC < 7 ) 

07500  PS1 (2,2>=APC(8> 

07600  UE(2,2)=APC(9> 

07700  WE(2»2)=APC( 10) 

07800  60  CONTINUE 

07900  70  CONTINUE 

08000  DO  75  1  =  1*  NXTL 

08100  75  READ ( TAPEVL  >  WNP 1 ( I ) *  WNPN ( I ) 

08200  IF( IFPRNT. EO.0)  GO  TO  100 

08300  REWIND  TAPEVL 

08400  WRITE <6. 230) 

08500  DO  90  K=1 *  NZTL 

08600  WRITE (TAPEOT, 240)  Z(K) 

08700  DO  e0  1=1, NXTL 

08800  READ (TAPEVL)  UEUF,WEUF 

08900  UEUF=UEUF~UFR5 

09000  WEUF»WEUF*UFRS 

09100  WRITE (TAPEOT , 250 )  I , X ( I ) , UEUF, WEUF 

09200  80  CONTINUE 

09300  90  CONTINUE 

09400  WRITE (TAPEOT, 252) 

09500  DO  95  1=1, NXTL 

09600  95  WRITE (TAPEOT, 250)  I , X ( I ) , WNP1 ( I ) , WNPN ( I ) 

09700  WR I TE ( TAPEOT , 255 ) 

09800  100  WRITE( TAPEOT, 180) 

09900  WRITE ( TAPEOT, 190)  TITLE 

10000  WR I TE ( TAPEOT , 200  >  NXSTRT, NZSTRT, NXT, NZT, NTR,  IFLOW,  ICHORD, 

10100  *  IFPRNT, ISPAN, ITMAX, IPZ* IPX, VGP,ETAE» 

10200  *  DETA( 1 ) , UFRS, PFRS, CNUFRS, EPS, EPSV, EPST 

10300  WRITE  (TAPEOT,  260)  SAL.FA,  SALFAS,  SBETA,  SBETAS,  SSIGMA,  SSGMAS, 

10400  *  SLAMDA, RW2, USTOP, Z IOTAE, Z IOTAL 

10500  WRITE( TAPEOT, 255 > 

10600  RETURN 

10700  120  FORMAT (9A6) 

10800  170  FORMAT ( 1P3E12. 4) 

10900  180  FORMAT (50X30HTHREE-D  BOUNDARY-LAYER  PROGRAM/ 1 H  , 

11000  *  47X36HFOR  TURBULENT  FLOW  ABOUT  A  SHIP  HULL) 

11100  190  FORMAT ( / / / /40X , 9A6 ) 

11200  200  FORMAT ( ///1H0, 32X0HNXSTRT=  , 13, 15X8HNZSTRT=  , 13, 15X8HNXT  =  ,13 

11300  *  / 1H0, 32X8HNZT  =  , 13, 15X8HNTR  =  , I 3 , 1 5 X8H I FLOW  =  ,13 

11400  *  / 1H0, 32X8HI CHORD=  , 13, 15X8HIFPRNT=  , 13, 15X8HISPAN  =  ,13 

11500  *  /1H0, 32 X8H ITMAX  =  ,I3,15X8HIPZ  =  ,I3,15X8HIPX  =  ,13 

11600  *  / 1H0, 32X7HVGP  =, 1  PE 1 4. 6 , 5 X7HETAE  =, E14. 6, 5X7HDETA  =,E14.6 

11700  *  / 1H0, 32X7HUFRS  =,£14.6, 5X7HPFRS  =, E14. 6, 5X7HNUFRS  =,E14.6 

11800  *  /1H0, 32X7HEPS  =, E14. 6, 5X7HEPSV  =, E14. 6, 5X7HEPST  =,E14.6> 

11900  210  FORMAT (IX, 25H**  INPUT  HULL  GEOMETRY  **//2X2HNZ , 1 X2HNX , 5X2HH1 . 

12000  *  10X2HH2, 10X2HK1, 10X2HK2, 10X3HK12, 9X3HK21 , 8X5HTHETA, 8X2HS1 / ) 
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12100  220  FORMAT ( 1H  t  2I3»  1P0E12. 4) 

12200  230  FORMAT ( 1 HI » 45H**  EXTERNAL  VELOCITY  DISTRIBUTION  (FT/SEC)  **/ 

12300  *  1H0,5X1HZ* 13X1HX» l3X2HUEt 12X2HWE/) 

12400  240  FORMAT ( 1H  »1PE12.5> 

12300  250  FORMAT ( 1H  * 7X» 13* 2X» 1P3E14. 5 ) 

12600  252  FORMAT ( // 1 H0i 19X 1HX i 12X4HWNP 1 » 1 1 X4HWNPN/ ) 

12700  255  FORMAT  < 1H1 ) 

12800  260  FORMAT < /////45X4SH**WILC0X-RUBESIN  TWO-EQUATION  TURBULENCE  MODEL* 

12900  *  ///45X29HTHE  CLOSURE  COEFFICIENTS  ARE : // 47X8H ALPHA  =,1PE12.3 

13000  *  5X8HALPHA*  =» E12. 3//47X8!  IBETA  =, E12. 3» 5X8HBETA*  =.E1 

13100  *  47X8HSIGMA  =» E12. 3. 5X8HSIGMA*  -?  E12. 3//47X8HLAMBDA  =»E1 

13200  *  5X8HRW2  =»E12.3// 

13300  *  /36X54HD I SS I  PAT I ON  RATE  PRESCRIBED  ANALYTICALLY  UP  TO  UPLUS 

13400  *  E12. 4///45X24HINITIAL  EDGE  CONDITIONS:// 

13500  *  47X8HI0TAE  E12..4, 5X8HI0TAL  =,E12.4) 

13600 


END 


ro  t  ) 


00100 
00200 
00300 
00400 
00500 
00600 
00700 
00800 
00900 
01000 
01  100 
01200 
01300 
01400 
01500 
01600 
01700 
01800 
01900 
02000 
02100 
02200 
02300 
02400 
02500 
02600 
02700 
02800 
02^00 
03000 
03100 
03200 
03300 
03400 
03500 
03600 
03700 
03800 
03900 
04000 
04100 
04200 
04300 
04400 
04500 
04600 
04700 
04800 
04900 


SUBROUTINE  LOG I C < NZ 1 , NZ2 ) 

INTEGER  TAPEIN* TAPEOT* TAPEGP* TAPEPE *  TAPEDT * TAPEVL 
COMMON/BLC0/  NX i NZ * NP i NXT* NZT , NTR , NPT , NXTL * N2TL *  NXSTRT » 

*  NZSTRT *  KC*  IT*  IFLOW,  I  CHORD »  ISPAN*  INTDIR 
C0MM0N/BLC9/  NSEP ( 41 ) * ICASE < 41 ) 

COMMON/TAPE/  TAPEINi TAPEOT » TAPEGP*  TAPEPF* TAPEDT * TAPEVL 

C  - - - - - - -  -  - - - - -  -  - 

NZ 1 P=NZ 1+1 
NZ2M=NZ2-1 

IF ( NX  *  NE , NXSTRT )  GO  TO  20 
DO  10  I =NZ 1 P i NZ2M 
10  I  CASE ( I ) “ 1 

I F ( I FLOW . EQ .  2 )  GO  TO  15 
I  CASE ( NZ 1 >=t 
I  CASE ( NZ2 ) =2 
RETURN 

15  I CASE ( NZ 1 ) *2 

I CASE ( NZ2 )  =  1 
RETURN 

20  I F  < I FLOW • EQ . 2 )  GO  TO  40 

I CASE ( NZ 1 ) =3 
I  CASE  <  NZ2 ) =5 
DO  30  I =NZ 1 P i NZ2M 
I  CASE < I )=4 

IF (NSEP ( 1  + 1  )  , NE. 0 )  I  CASE ( I ) ^6 
IF ( NSEP < 1-1 ) .NE.0)  I  CASE < I > =3 

IF ( NSEP  < 1+ 1 ) . NE . 0. AND . NSEP ( I  —  1 > . NE -  0 )  GO  TO  25 
GO  TO  30 

2*5  NSEP  <  I  )  .  NF .  0  >  GO  TO 

NSEP  < I >  =  1 

WRITE ( TAPEOT t  100)  I  *  NX 
30  CONTINUE 
RETURN 

40  I  CASE ( NZ2 ) =3 
I  CASE  < NZ 1 ) =5 
DO  50  I =NZ 1 P *  NZ2M 
I  CASE ( I ) «4 

IF ( NSEP < 1-1  >  . NE, 0)  I  CASE ( I ) =6 
IF ( NSEP < 1+1 ) .NE.0)  I  CASE < I ) =3 

IF  <  NSEP  < I  + 1 ) • NE, 0- AND • NSEP  < I- 1 ) • NE, 0 )  GO  TO  45 
GO  TO  50 

45  IF < NSEP < I ) , NE, 0)  GO  TO  50 
NSEP ( I )=1 

WR I TE  <  TAPEOT  *  100)  IiNX 
50  CONTINUE 
RETURN 

100  FORMAT  < /9X26H**  CALCULATIONS  ALONG  NZ  14, 2X 1 8HTERMINATED  AT  NX 

*  , I4/9X40H  BECAUSE  OF  MERGING  SEPARATION  REGIONS//) 

END 


00100 

00200 

00300 

00400 

00500 

00600 

00700 

00800 

00900 

01000 

01100 

01200 

01300 

01400 

01500 

01600 

01700 

01800 

01900 

02000 

02100 

02200 

02300 

02400 

02500 

02600 

02700 

02800 

82900 

03000 

03100 

03200 

03300 

03400 

03500 

03600 

03700 

03800 

03900 

04000 

04100 

04200 

04300 

04400 

04500 

04600 

04700 

04800 

04900 

05000 

05100 

05200 

05300 

05400 

05500 

05600 

05700 

05800 

05900 

06000 


20 


80 


COMMON/BLC3/ 
COMMON/BLC5 / 


SUBROUTINE  OUTPUT (LL> 

DIMENSION  Y< 101 ) *  SUMTBL ( 10*41)* CVC( 101 ) *  BETA ( 101 ) * USUSE< 101 ) 
DIMENSION  WPRNT  < 101 ) * ELPRNT ( 101  > *SVC( 101 ) 

INTEGER  TAPE IN* TAPEOT*  TAPEGP* TAP EPF* TAPEDT * TAPEVL 
COMMON/ BLC0/  NX*  NZ* NP* NXT  * NZT  *  NTR* NPT  *  NXTL  *  NZTL  *  NXSTRT  * 

NZSTRT *  KC *  IT*  IFLOW*  ICHORD*  ISPAN*  INTDIR 
X<81 ) i Z<41 ) *  DETA ( 101 ) *  ETA ( 101 ) *  ETAE *  VGP* CEL *  BEL 1* BEL 
F( 101*2*2) *U( 101*2*2) *V( 101 )i6( 101*2*2) *W( 101*2*2)* 

T ( 101 ) v  TPROF ( 620 ) *  TPCF ( 10) 

NSEP(h 1  ) *  I  CASE (41  > 

PS1 <2*  2) *  HI <2*  2) *  H2C2*  2) *  UE(2*  2) * 

WE (2*2) *  CK 1(2*2)* CK2<2* 2 ) *  CK 1 2 ( 2 *  2 ) * CK21 (2*  2) * 

THETA ( 2*  2 ) * PFRS*UFRS*  CNUFRS * UREF 1  * WNP 
EDV( 101 ) 

SALT  A  * SALFAS*  SBETA*  SBETAS*  SLAMDA*  SSIGMA*  SSGMAS*  USTOP- 
Z I OT AE »  Z I OTAL  *  RW2 
CON*  EPS*  ICONVE*  ITMAX*  WTEDG 

E( 101  *  2*  2) * WT< 101 ) * WT2 ( 101*2*2)*  ELT ( 101 ) 

DUDX  *  DWDX  *  DUDZ*  DWDZ 
IPZ*  IPX*  I  PR I NT  *  EPSV*  EPST 

TAPEIN*  TAPEOT, TAPEGP * TAPEPF *  TAPEDT *  TAPEVL 


COMMON /BLC9/ 
COMMON/ PRSC/ 


COMMON /BLCE/ 
COMMON/CONS/ 


COMMON/STOR/ 
COMMON/TURe/ 
COMMON/DERV/ 
COMMON/ I PRT/ 
COMMON/TAPE/ 
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:)* 


* ( UE ( 2 *  2) /HI <2* 2 
2 ) +UE (2*2)* 

*DWDX+WE(2*  2) /H2(2*  2)* 
*  2  )  )  **  l  Ub  l  2  *  )  * 

*2) ) ) /SIN ( THETA ( 

2(2*2) *  COS ( THET A 


:)  ) 


GO  TO  (10*220)*  LL 

USE=SORT ( UE  <  2* 2 > **2+WE ( 2* 2 > **2+2. *UE ( 2* 2 ) *WE (2*2)* 

*  COS (THETA (2*  2) ) ) 

DUSEDS= ( ( UE (2*2) +WE (2*2) *COS ( THET A ( 2 »  2  >  )  ) 

*  DUDX+WE ( 2 *  2) /H2(2*  2)*DUDZ >+(WE(2* 

*  COS (THETA (2*  2)  )  > * (UE ( 2*  2 ) /HI (2*  2) 

*  DWDZ  )  -»-UE  <  2*  2  ;  *wt  <  2*  ;  *6i  N  k  i  Ht  I  a  ( 2 

*  ( <  CK 12(2*2)  +  CK2 1 ( 2* 2 ) *COS ( THETA ( 2 

*  CK1 (2*2)  >+WE(2*2>*( (CK21 (  2*  2M  CK1 

*  SIN (THETA (2* 2) >+CK2(2*2> ) ) ) /US E** 

URUE*UREF 1 /UE (2*2) 

US  1  =  (  UE  (  2*  2  )  /USE  )  **2 
US2=  (  UREF  1  /USE  )  *-*2 

US3=2 *  *UE (2*2) *UREF 1 /USE**  2*C0S ( THETA (2*2 
B 1 =URUE*S IN ( THET A ( 2 *  2 )  ) 

B2«SURUE*  COS  (  THET  A  ( 2  *  2 )  ) 

RX=UE ( 2 *  2 ) *PS 1 (2*2) /CNUFRS 
TRCMN=SQRT( CNUFRS*PS1 ( 2* 2 ) /UE ( 2* 2 ) ) 

WS I NTH=B 1 /URUE 
XFCKC.EQ. 1 )  WS I NTH=0 . 

NPYES=0 

NPM03=(NP-1 ) /3 
ENPM03=(NP-1 . ) /3. 

D I FFNP»ENPM03— NPM03 
IF ( DIFFNP. GT. . 1 )  NPYES=1 
DO  20  J= 1  *  NP 
Y ( J ) =TRCMN*ETA ( J ) 

CVCT =UREF  1  /  UE  (  2 *  2  )  *(.  S(THETA(2*2)  ) 

CVCD«( 1 • +WE (2*2) /UE (2*2) *COS ( THETA ( 2*2) ) ) *URUE 
DO  80  J= 1  *  NP 

USUSE ( J ) —SORT ( US 1 *U ( J  *  2  *  2 ) **2+ US2*W ( J  *  2  *  2 ) **2+US3*U ( J  *  2 *  2 ) * 

*  W  (  J  *  2  ?  2  ) ) 


2*2) >  + 
(2*2)))/ 


BETA ( J ) =0. 

IF ( J • GT • 1 . AND. KC. NE. 1 ) 

*  BETA ( J ) =57 . 295 78* AT  AN ( B 1 *W (J*2*2)/(U(J*2*2) +  B2*W ( J  *  2*2))) 
SVC ( J ) =WS I NTH*W ( J  *  2*  2 ) 

CVC( J)«(U( J*  2*  2)+W( J* 2*2) *CVCT ) / CVCD 

IF ( KC • NE . 1 )  BETA ( 1 ) =57, 295 78* A TAN ( B 1 *T ( 1 )/(V( 1 )+B2*T( 1 ) ) ) 


-76- 


0610® 

06200 

06300 

06400 

06500 

06600 

06700 

06800 
06900 
07000 
07100 
07200 
07300 
07400 
07500 
07600 
07700 
07800 
0790© 
08000 
08100 
08200 
08300 
08400 
08500 
08600 
08700 
08800 
06900 
09000 
09100 
09200 
09300 
09400 
09500 
09600 
09700 
09800 
09900 
10000 
10100 
10200 
10300 
10400 
10500 
10600 
10700 
10800 
10900 
1  1000 
11100 
11200 
11300 
11400 
11500 
11600 
11700 
11800 
11900 
12000 


IF  < I PRINT . EO,  0 )  GO  TO  111 
WRITE <  TAPEOT *  280 ) 

WSCALE=US1 / ( US2*RX ) 

ELSCAL=URUE*PS1 (2*2) 

WPRNT ( 1 ) =9. 99999E+37 
ELPRNT ( 1 )=0. 

DO  84  J=2>NP 
WPRNT ( J ) =WT ( J ) *WSCALE 
84  ELPRNT  <  J  )  =EL.T  (  J  )  *ELSCAL 

WRITE (TAPEOT, 290)  (JiY(J), CVC ( J) i SVC< J> *  USUSE ( J > i BETA ( J ) * 

*  E( Ji2*2> i WPRNT ( J) i ELPRNT ( J) , EDV ( J ) *  J=l»NPi3) 

IF ( NPYES. EO. 0 )  GO  TO  111 

WRITE < TAPEOT i 290)  NP. Y(NP) * CVC(NP) , SVC(NP) > USUSE (NP) i BETA ( NP  > 

*  E(NP,2*2) * WPRNT (NP) * ELPRNT (NP) iEDV(NP) 

1 1 1  BETA < NP ) =BETA ( NP ) /57. 29578 

CTRM=2. /SORT ( RX ) 

UEUF=UE ( 2 i 2) /UFRS 
S 1 SQRX=PS 1 (2i 2) /SORT ( RX ) 

COSTH=COS( THETA (2i 2) ) 

CTRM2=UE ( 2 i 2 ) *UREF 1 /UFRS**2*T ( 1 ) *COSTH 
CFC=CTRN*  (  UEU F**2*V  (  1  )  +CTRM2 ) 

CFN=CTRM*UE(2i2)*UREFl/UFRS**2*T< 1 ) *SIN < THETA ( 2, 2) ) 

C I 2=0 . 

C I 3=0 . 

Cl 4=0 . 

WEUECT =WE ( 2 1 2 ) /UE (2,2) *COSTH 
URUECT =UREF 1 /UE ( 2 , 2 ) *COSTH 
C2=(U( 1 i 2 1 2)+URUECT*W( 1 , 2 *  2 )  ) **2 
C3=w < i ,  2,  2  >**2 
C4=U ( If 2»2)**2 
DO  170  J=2iNP 

C22=(U(  Ji  2,  2)+URUECT*W< J,  2»  2)  )  **2 
C33=W< Ji 2f 2)**2 
C44=U< Jf 2f 2>**2 
C 1 2  =  C 1 2  +•  • 5*  <  C2+C22 ) *DETA  < J— 1 ) 

CI3=CI3+. 5*< C3+C33 )*DETA( J~1 ) 

CI4=CI4+, 5*( C4+C44 )*DETA( J-l  ) 

C2=C22 
C3=C33 
170  C4=C44 

TNUM=F(NPi 2* 2)+URUECT*G(NP, 2i 2) 

TDEN= 1 • +WEUECT 

DLSTS=S 1 SORX  *  <  ETA  <  NP )  —  ( F (NP» 2, 2) +UREF 1 /UE ( 2» 2)*G(NP, 2, 2) 

*  *COSTH) /TDEN) 

THTAS=S1S0RX* ( TNUN/ TDEN—  C I 2/TDEN**2 ) 

IF ( KC . GT, 1 )  GO  TO  180 
CFC=CTRM*UEUF**2*V(  1  ) 

CFN=0 . 

DLSTS=S 1 SORX* (ETA(NP)-F (NP* 2 ,2  >  > 

THTAS=S 1SORX* (F(NPi 2i 2 )-CT4 ) 

180  IF< ABS(WE(2, 2) ) .  GT. 1 . E-8 )  GO  TO  200 
DLSTN=0. 

THTAN=0. 

GO  TO  210 

200  DLSTN=S1 SORX *( ETA (NP)-UREFl /WE(2i 2>*G(NPi 2i 2) ) 

THT AN  =  S 1 SORX  *UREF 1 /WE ( 2 f  2 ) * ( G ( NP »  2 1 2  > -UREF 1 /WE  <  2f 2 ) *C 1 3 ) 

210  IF (  I  PRINT, EQ.  1  )  WR I TE ( TAPEOT i 300 )  CFC t  DLSTS »  THTAS » BETA ( 1  ) i CFN 

*  DLSTNt THT AN » RXiUSEi DUSEDS 
SUMTBL ( 1 i NZ ) — V ( 1 ) 

SUMTBL (2*NZ)=T( 1 ) 


12100 

12200 

12300 

12400 

12500 

12600 

12700 


SUMTBL ( 3, NZ ) —  CFC 
SUMTBL  <  4 1 NZ ) =  CFN 
SUMTBL  <  5 i NZ ) =BET A < 1 ) 
SUMTBL (6iNZ) =DLSTS 
SUMTBL  <  7,  NZ  )--DLSTN 
SUMTBL  <  8, NZ ) =THTAS 
SUMTBL  <  9i NZ ) -THTAN 


12800 

12900 

13000 
13100 
13200 
13300 
13400 
13500 
13600 
13700 
13800 
13900 
14000 
14100 
14200 
!  4300 
14400 
14500 
14600 
14700 
14800 
1  4900 
15000 
15100 
15200 
15300 
15400 


SUMTBL ( 10, NZ )=RX 
RETURN 

220  WRITE  <  TAPEOT ,  310  >  X(NX> 

DO  240  I=NZSTRT , NZT 

IF  <  NSEP (  I  )  • EO. 0 )  GO  TO  230 

WRITE (TAPEOT, 340)  I, Z( I ) 

GO  TO  240 

230  WRITE ( TAPEOT, 320 )  I , Z ( I ) , ( SUMTBL (J, I ) , J=1 * 10) 

240  CONTINUE 

WRITE (TAPEOT, 330) 

RETURN 

C  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  ~  —  —  —  —  —  —  — 

280  FORMAT ( 1 H0 , 2X1 HJ ,4X1 HY , 1 0X4MU/UR, 9X4HW/UR , 9X4HO/OE , 9X4HBETA , 

*  8X7HE/UR**2, 6X5H0MEGA , 9X2HEL , 1 1X3HEPS) 

290  FORMAT ( 1H  ,  13*1  PE 10. 3, 8E 13.4) 

300  FORMAT (2X6HCFC  1  PEI  4. 6, 3X6HDLSTS=, El  4. 6, 3X6HTHTAS-, El  4. 6, 

*  3X6HBET  A 1  = ,  E 1  4 . 6/2X6HCFN  =, El 4. 6, 3X6HDLSTN=, E 1 4, 6 , 

*  3X6HTHTAN=, E14. 6, 3X6HRX  - , El 4 . 6/2X6H0E  =, E14- 6, 

*  3X6HD0EDX=  ,  E 1  4 . 6/1H0, 62  (  2H*-—  )  ) 

310  FORMAT  (  1H0,  35  X28H *•#-*•* —  SUMMARY  TABLE  FOR  X  =*1PE10.2, 

*  1 X7H  -- .  »***/iM0* 3H  NZ , 5  X 1 HZ , 

0X5 HV WALL •,  6X5HTNAL.L  ?  ^XZHPFC,  0X3HCFN-  7y*,L,PETA* 

*  7X5HDLSTS, 6X5HDLSTN, 6X5HTHTAS, 6X5HTHTAN, 7X2HRX ) 

320  FORMAT ( 1H  *13,  1  PEI 0.3*  10E11.3) 

330  FORMAT ( 1 H0 , 62 ( 2H*- ) / ) 

340  FORMAT  (  1  H  ,13,1  PE  10.3,3X1 7H*-**  SEPARATED  ***> 

END 


00100 
00200 
00300 
00400 
00500 
00600 
00700 
00300 
00900 
01000 
01100 
01200 
01300 
01400 
01500 
01600 
01700 
01800  C  - 
01900 
02000 
02100 
02200 
02300  10 

02400 
02500  20 

02600  30 

02700 
02800 
02^00 
03000 
03100 
03200 
03300 
03400 
03500 
03600 
03700 
03800 
03900 
04000 
04100 
042:00 
04300 
04400 
04500 
04600 
04700 
04800 
04900 
05000 
05100 
05200 
05300 
05400  40 

05500 
05600 
05700 
05800 
05900 
06000 


SUBROUTINE  PRESUR 

INTEGER  TAPE IN* TAPEOT » TAPEGP*  TAPEPF* TAPEDT * TAPEVL 
DIMENSION  ARC (10) 

COMMON /BLC0/  NX  *  NZ  *  NP,  NXT*  NZT*  NTR*  NPT  *  NXTL. *  NZTL  *  NXSTRT  * 

*  NZSTRT iKC» IT* IFLOW* ICHORDi ISPAN* INTDIR 
C0MM0N/BLC3/  X  <81 >  »  Z  <41 ) *  DETA  < 101 ) *  ET A ( 101 ) *  ETAE»  VGPi CEL i BEL1 *  BEL 
COMMON/ CHRD/  WNP 1<81)*  WNPN  <  8 1  ) 

COMMON/ PRSC/  PS l(2i2)iHl(2i2)i H2  <2i2)i UE (2»2)i 

*  WE <2*  2) i CK1 <2*  2) * CK2C2* 2) i CK12<2«2) * CK21 <2*  2) * 

*  THETA  <  2*2) i PFRStUFRS*  CNUFRS * UREF 1  * WNP 
COMMON/ PRES/  P 1 < 2*  2 >  »  P2 < 2 *  2 ) i P3 < 2i 2 ) *  P4 ( 2*  2 ) *  P5 ( 2*  2 ) i 

*  P6<  2*  2) *  P7(2*  2) *  PB(2f 2) *  P9(2*  2) *  P10(2*  2) *  PI  1 <2*  2) * 

*  P12<  2*  2) *  P13< 2*  2) *  P14(2v2) 

COMMON/DERV/  DUDX* DWDX*  DUDZ*  DWDZ 

COMMON/TAPE/  TAPE  IN*  TAPEOT*  TAPEGP*  TAP EPF* TAPEDT*  TAPEVL 
COMMON/ I PRT /  IPZ*  IPX*  I  PR I NT , EPSV * EPST 
COMMON /SVNM/  IPCF* IPRF 


IF(NZ.EO.NZT)  GO  TO  20 
IF(NZ. 6T. NZSTRT )  GO  TO  10 
ICNT=(NZ~1 )*NXTL+NX 
GO  TO  30 

I CNT  = ( NZ-2 ) *NXTL+NX 
GO  TO  30 

I CNT = ( NZ~3 ) *NX  TL+NX 
I PNT2= I CNT 

CALL  STORIT (TAPEGP* APC* IPCF* IPNT2* 1 ) 
HI  1 =APC  < 1 ) 

H:?t=APC<2> 

THET A 1 -  APC  <  7 ) 

SI  1 =APC ( 8 ) 

UE 1 -APC ( 9 ) 

WE 1 =APC  <  10) 

ICNT=ICNT+NXTL 
I  PNT*2=  I  CNT 

CALL  STORIT  <  TAPEGP i APC i  IPCF*  IPNT2*  1  ) 
H12=APC( 1 ) 

H22«APC<2) 

THETA2=APC(7) 

S12=APC(B) 

UE2=APC  <  9 ) 

WE2=APC< 10) 

ICNT=ICNT+NXTL 
I  PNT2==  I  CNT 

CALL  STORIT (TAPEGP* APC* IPCF* IPNT2i 1 ) 
H 1 3=APC ( 1 ) 

H23=APC<  2) 

THET A3=APC ( 7 ) 

S 1 3=APC ( 3 ) 

UE3“APC<9) 

WE3— APC  <  10) 

IF(NZ. EO. NZT)  GO  TO  60 
IF  <  NZ  • GT. NZSTRT )  GO  TO  50 
Z  1  =*Z  (  NZ  ) 

Z2=Z  <NZ+1  ) 

Z3=Z  <NZ  +  2) 

A1«<Z1-Z2>*( Z1-Z3) 

A2=( Z2-Z1 )*( Z2-Z3) 

A3= ( Z3— Z 1 ) * ( Z3- Z2 ) 

D1  *=  (  2.  *Z1-Z3-Z2> /A1 


06100 

D2=(Z1~Z3)/A2 

06200 

D3“ ( Z 1 — Z2 ) / A3 

06300 

DWDZ=UREF1*WNP1 (NX) 

06400 

GO  TO  70 

06500 

50 

Z 1 —Z ( NZ—  1 ) 

06600 

Z2--Z  <  NZ  ) 

06700 

Z3=Z  <  NZ  + 1 ) 

06800 

A 1  =*  (  Z2-Z  1  )  *  <  Z3—Z  1  ) 

06900 

A2=<  Z2-Z1 )*<Z3-Z2> 

07000 

A3=(Z3-Z2)*(Z3-Z1 ) 

07100 

D1=-(Z3-Z2)/A1 

07200 

D2- ( Z3-2. *Z2+Z1 ) / A2 

07300 

D3=<  Z2-Z1 ) / A3 

07400 

DWDZ=D1*WE1+D2*WE2+D3*WE3 

07500 

GO  TO  70 

07600 

60 

Z1=Z (NZ-2) 

07700 

Z2=Z ( NZ— 1 ) 

07800 

Z3“Z ( NZ ) 

07900 

A1=(Z2-Z1 )•*(  Z3-Z1 ) 

08000 

A2= ( Z2— Z 1 ) * ( Z3—Z2 ) 

08100 

A3=<  Z3-Z2 >* ( Z3-Z1 ) 

08200 

Dl  =  (  Z3-Z2) /A1 

08300 

D2=—  <  Z3—Z 1 ) / A2 

08400 

D3=  <2- *Z3-Z1-Z2)/A3 

08500 

DWDZ=UREF 1 *WNPN ( NX ) 

08600 

70 

DUDZ=D1*UE1+D2*UE2+D3*UE3 

08700 

WNP=DWDZ/UREF1 

08800 

DH1  DZ  =  D  1  *H  1  1  +  D2-*H  1  2-*  D3*H  1 3 

08900 

DTmjj  Z  =D  i  *THET a  i  +  D2*-  i  ml  i  i  Mt.  t  ftv 5 

09000 

DS1 DZ~ Dl *S1 1 +D2*S 1 2+D3*S 1 3 

09100 

CPR1~SQRT( CNUFRS*UE1*S1 1 ) *H21*SIN ( THETA 1 >*H1  1 /H21 

09200 

*  *UREF1/UE1 

09300 

CPR2— SORT ( CNUFRS*UE2*S 1 2 ) *H22*S I N ( THETA2 ) *H 1 2/H22 

09400 

*  *UREF1/UE2 

09500 

CPR3=S0RT ( C,NUFRS*UE3*S 13) *H23*SIN ( THETA3  >  *H13/H23 

09600 

*  *UREF1/UE3 

09700 

DCPDZ=D1*CPR1+D2*CPR2+D3*CPR3 

09800 

IF (NX. EQ. NXT)  GO  TO  120 

09900 

IF ( NX . GT . NXSTRT  >  GO  TO  110 

10000 

100 

ICNT=(NZ-1 ) *  NXTL+NX 

10100 

I  NX  =  1 

10200 

GO  TO  130 

10300 

110 

I CNT — ( NZ— 1 ) *NXTL+NX— 1 

10400 

INX=2 

10500 

GO  TO  130 

10600 

120 

I CNT  = ( NZ— 1 )*NXTL+NX-2 

10700 

I  NX=3 

10800 

130 

I PNT2= I CNT 

10900 

CALL  STOR I T ( TAPEGP i APC* IPCF» IPNT2i 1 ) 

1  1000 

H2 1 =APC ( 2 ) 

1  1100 

THETA 1 =APC ( 7 ) 

1  1200 

SI  1 =APC ( 8 ) 

1  1300 

UE 1 =APC ( 9 ) 

11400 

WE 1 =APC ( 10) 

11500 

I CNT  « I CNT  + 1 

11600 

I PNT2= I CNT 

11700 

CALL  STOP I T ( TAPEGP »  APC  i  IPCFi  IPNT2i  1 ) 

11800 

H22=APC(2) 

11900 

THETA2=APC ( 7 ) 

12000 

S 1 2— APC ( 8 ) 

12100 

12200 

12300 

12400 

12500 

12600 

12700 

12800 

12900 

13000 

13100 

13200 

13300 

13400 

13500 

13600 

13700 

13800 

13900 

14000 

14100 

14200 

14300 

14400 

14500 

14600 

14700 

14800 

14900 

15000 

15100 

15200 

15300 

15400 

15500 

15600 

15700 

15800 

15900 

16000 

16100 

16200 

16300 

16400 

16500 

16600 

16700 

16800 

16900 

17000 

17100 

17200 

17300 

17400 

17500 

17600 

17700 

17800 

17900 

18000 


UE2=APC  <  9 ) 

WE2=APC  < 10) 

I CNT-I CNT+ 1 
I PNT2= I CNT 

CALL  STORIT  <  TAPEGP*  A PC* IPCF* IPNT2* 1 ) 

H23==APC  <  2 ) 

THETA3=APC ( 7 ) 

S13=APC(8) 

UE3—APC ( 9 ) 

WE3—APC ( 10) 

IF ( INX • EO- 3 )  GO  TO  160 
IF < INX . EQ. 2 )  GO  TO  150 
140  X 1  =  X  ( NX ) 

X2=X(NX+1) 

X3=X (NX+2) 

A1=(X1-X2)*(X1-X3) 

A2=( X2-X1 )*(  X2-X3) 

A3=  <  X3— X 1 )*(X3-X2) 

Dl  =  (2. *  X 1— X3—X2 ) / A1 
D2=(X1-X3)/A2 
D3=(  X1-X2) /A3 
GO  TO  170 
150  X 1 =X ( NX— 1 ) 

X2-X (NX ) 

X3=X  (  NX  + 1  ) 

A 1  =  ( X2—X 1 >*(X3-X1 ) 

A2= ( X2— X 1 )*(X3-X2> 

A3= ( X3-X2 ) *  <  X3— X 1 ) 

Di =— <  X3-X2 > / A i 
D2=(X3-2.*X2+X1 )/A2 
D3=( X2-X1 ) /A3 
GO  TO  170 
160  X  1  =  X  < NX— 2 ) 

X2=X (NX-1 ) 

X3=X ( NX ) 

A1  =  <  X2— X  1  )  -*  (  X3--X  1  ) 

A2=( X2-X1  )*<  X3-X2) 

A3=  (  X3— X2  )  *  (  X3— X  1  ) 

D1  =  (  X3-X2 ) / A 1 
D2=— ( X3-X 1 ) /A2 
D3= ( 2 . ■*  X3— X 1  —  X2  ) / A3 
1 70  DUDX=D 1 *UE 1  +  D2-*UE2+D3*UE3 

IF(NX.EQ.NXSTRT)  DUDX- ( UE2-UE 1 )/(X2-Xl ) 
DBPDX=D1*SQRT< CNUFRS*UE 1 *S 1 1  > *H21*S IN ( THETA 1 )  + 

*  02* SORT ( CNUFRS*L'E2*S 1 2 ) *H22*S I N ( THET A2 ) + 

*  D3*SQRT ( CNUFRS*UE3*S 1 3 ) *H23*S I N ( THET A3 ) 

1 90  DWD  X  =D 1 *WE 1 +D2*WE2+D3*WE3 

HTRM1=H21 *S IN ( THETA 1 ) 

HTRM2- H22-*S I N ( THETA2 ) 

HTRM3=H23*S IN ( THET A3 ) 
DHTRMX=D1*HTRM1+D2*HTRM2+D3*HTRM3 
COTTH=COTAN<  THETA (2. 2)  ) 

SINTH=S IN  <  THETA (2*2)  ) 

COS TH— COS  <  THETA (2* 2 )  ) 

BPR^SQRT ( CNUFRS  *UE (2*2) *PS1 (2*2)  ) *H2 (2*2) *  S I NTH 
H1H2=H1 (2*2) /H2 (2*2) 

IF ( KC  * EQ . 4 )  GO  TO  460 
DUDZ=0. 

DWDZ=0. 

DS1DZ=0. 

-81- 


18100  DCPDZ=0. 

18200  460  DL0GS1 =DS 1 DZ /PS 1 ( 2»  2 ) 

18300  PI  (  2»  2  )  =. 5*( PS1 (2.2) / (HI (2i 2>*UE(2i 2) )*DUDX+1 .  )+PSl <2i 2) / 

18400  *  (HI (2. 2 ) *H2<2»  2)*SINTH ) *DHTRMX 

18500  P2 (2*2) =PS1 <2*2>/<UE(2*2>*Hl (2.2) >*DUDX-PS1 <2*2>*CK1 (2.2>*C0TTH 

18600  P3(2.2)=-PS1 (2.  2)*C0TTH*CK2(2. 2 > *UREF1 /UE ( 2. 2) 

18700  P4 ( 2t  2 ) =PS 1 (2*2) *CK21 (2.2) 

18800  P5(2»2)=PS1 <2.2)/H2<2*2)*UREFl/UE<2.2>**2*DUDZ+CK12<2. 2)* 

18900  *  PS1 (2»2>*UREF1/UE(2.2> 

19000  P6 (2.2) =PS 1 (2. 2) /(HI ( 2* 2 ) *BPR > *DCPDZ 

19100  P7(2.2>~PS1 (2»2)/H2(2*2)*UREFl/UE<2»2> 

19200  P8 (2.2) =PS 1 ( 2»  2 ) *CK2 ( 2t  2 ) /S INTH* ( UREF 1 /UE (2.2)  )**2 

19300  P9 ( 2 1 2 ) =PS 1 (2*2)* CK 1 (2.2) /SINTH*UE (2*2) /UREF1 

19400  P10(2*  2)=PS1 (2* 2) /HI (2*2) 

19500  Pll (2.2) =PS1 <2»2>*( 1. /(UE(2*2)*H1 (2*2) > *DUDX+WE<2* 2 > / (UE <2* 2 > 

19600  *  **2*H2(2.  2) ) *DUDZ-COTTH*CKl ( 2. 2)+CK2(2*  2) /SI NTH* 

19700  *  (WE(2*2)/UE<2.2> >**2+CK12(2. 2>*WE(2. 2) /UE(2* 2) ) 

19800  P 1 2 (  2. 2  > =PS 1 ( 2*  2 ) / ( UE (2*2) *UREF 1 > * ( UE ( 2* 2 > /HI (2.2>*DWDX+ 

19900  *  WE( 2*  2) /H2( 2*  2>*DWDZ— C0TTH*CK2 (2*2) *WE (2*  2 )**2+ 

20000  *  CK 1(2*2) /SINTH*UE ( 2 *  2 ) **2+CK21 (2*  2 ) *WE ( 2* 2 ) *UE ( 2* 2 ) ) 

20100  P13(2»  2>=1 . -PS1 (2*  2) / (HI (2*  2)*UE(2*  2) >*DUDX 

20200  PI 4 <2*  2) “P7 ( 2*  2)*( DLOGS1— DUDZ/UE(2* 2) ) 

20300  IF(KC.NE.l)  RETURN 

20400  DWDZ=UREF 1 *WNP 

20500  NX1=NX+1-INX 

20600  NX2=NX 1+1 

20700  NX3=NX2+1 

20800  IF  ( NZ  .  E<5.  NZSTRT  )  D2WDZX=D1*WNP1  (NX1  )+D2*WNPl  ( NX2 ) +D3 *WNP1  ( N  X  3  ) 

2 0900  IF<NZ.  EO.  NZT)  B2NDZX=*D1*'.*,N°N(NX1  )  +D2*WNPN  (  NX2 '  +n,3*(JNP!\(  ( NX3 ) 

21000  P3 (2*2) ~P7 (2*2) 

21100  P5 (2*2) =0 . 

21200  P6(2. 2>=P7(2. 2) 

21300  P8 (2*2) =0 . 

21400  P9 ( 2i 2  >  =0. 

21500  P12 ( 2*  2 ) =PS1 (2.2) /UREF 1 * ( UREF 1 *D2WDZ  X /HI (2*2) 

21600  *  +DWDZ**2/(UE(2»2)*H2(2»2) >+CK21 (2* 2>*DWDZ > 

21700  RETURN 

21800  END 
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SUBROUTINE  PROFIL(LL) 

COMMON /BLC0/  NX»NZ»NPi NXT, NZT, NTR, NPT ,  NXTL , NZTL, NXSTRT, 
*  NZSTRT ,  KC ,  IT,  IFLOW,  I  CHORD,  ISPAN, INTDIR 


COMMON/ BLC3/ 
COMMON /BLC5/ 

* 

COMMON /BLCE/ 
COMMON/ PRSC/ 

* 

* 


X(81)vZ(41 ) , DETA  < 101 ) , ETA  < 101 > , ETAE , VGP , CEL , BEL 1 , BEL 
F(101,2,2),U(101,2,2) , V( 101 ) ,  G< 101,2,2) ,W( 101,2,2), 
T( 101 ) , TPROF ( 620 ) , TPCF ( 10) 

EDV  < 101 ) 

PS1 (2,2) ,H1 (2,2) ,H2(2t 2) ,UE(2, 2) , 

WE (2,2) , CK1 (2,2) , CK2(2, 2) , CK12(2, 2)  ,  CK21 (2,2)  , 

THETA ( 2,2) , PF RSiUFRS, CNUFRS,UREF1 , WNP 


COMMON/TURB/ 
COMMON/STOR/ 
COMMON/ CONS/ 


Ed  01, 2,  2),  WT  (101),  WT2<  101,2,2),  ELT  (  101  > 

CON, EPS,  ICONVE,  I  TMAX  ,  WTEDG 

SALFA, SALFAS, SBETA, SBETAS, SLAMDA, SSIGMA, SSGMAS, USTOP 
ZIOTAE,  ZIOTA,  ,  RW2 


GO  TO  (10,40,70),  LL 

IF (NTR, EQ. 0 )  GO  TO  15 

CALL  START 

GO  TO  40 

V( 1 >=.332 

T ( 1 )=. 332 

WCON=20 . /SBETA 

REY=UE(2, 2)*PS1 (2, 2) /CNUFRS 

WTEDG=UREF1/UE<2, 2 > *SQRT ( REY*Z I OTAE/SALFAS > / ( Z IOTAL*ETAE ) 
C1=ETA(NP)*V( 1 ) 

C2=3. -2. *C1 
C3=— 2. +  C1 
ClD2=Cl/2. 

C3D4=C3/4. 

DO  20  J=  1 , NP 

ERAT =ET A ( J  > /ETA ( NP ) 

ERAT  C2=ERAT*C2 
ERATSQ-ERAT  *  *2 
ERTSQ3=ERATSQ*C3 

F  <  J ,  2,  2  )  *  (  C 1 D2+ERATC2/3.  +C3D4~ERATS0 )  *ERATSO*ET A  (  NP  ) 

U( J, 2, 2)  =  ( C1+ERATC2+ERTSQ3)*ERAT 
G(J,2,2)=F( J,2,2) 

W( J,2i 2)=U< J, 2,2) 

E( J, 2, 2>«0. 

IF(J.EQ.l)  GO  TO  20 

WT ( J ) -WCON/ETA ( J ) **2 

IF ( WT ( J ) . LT. WTEDG)  WT( J)=WTEDG 

WT2 ( J , 2  v  2 ) =WT ( J ) **2 

CONTINUE 

E(NP,2,2)=ZI OT AE 
WT2 (1,2,2) =WT2 (2,2,2) 

K=NP+ 1 
L=NPT 

DO  60  J =K , L 

F  <  J ,  2 , 2 ) =ET A  <  J  >  +F  <  K- 1 , 2 , 2  > -ET A  <  K- 1 > 

U  (  J ,  2 , 2  )  =  1  . 

W< J,2,2)=W(K-1,2,2) 

G  ( J,  2,  2  )  =W  (  K—  1 , 2, 2 )  *  (  ETA  <  J  )  —  ETA  <  NP  >  )  +G  (  K— 1,2,2) 

E  <  J , 2 , 2 ) =E ( K— 1,2,2) 

WT ( J ) =WTEDG 
WT2( J, 2, 2 ) =WTEDG**2 
IF ( LL, NE . 1 )  RETURN 
CALL  ARRYMV (4, 1,2) 

RETURN 
K=NP+1 


00100  SUBROUTINE  SHIFT < NZ 1 , NZ2 > 

00200  INTEGER  TAPEIN, TAPEOT , TAPEGP, TAPEPF ,  TAPEDT ,  TAPEVL 

00300  COMMON/BLC0/  NXiNZiNPiNXTiNZTi NTR,  NPT i NXTLi NZTL , NXSTRT , 

00400  *  NZSTRT, KC, IT, IFLOW, ICHORD, ISPAN, INTDIR 

00500  C0MM0N/BLC5/  F(101,2,2),U(101,2,2>,V(101>,G(101»2,2),W(101,2,2>, 

00600  *  T( 101 > , TPROF(620> , TPCF( 10) 

00700  C0MM0N/BLC9/  NSEP < 4 1 ) ,  I  CASE ( 41 ) 

00000  COMMON/TAPE/  TAPEIN, TAPEOT, TAPEGP* TAPEPF, TAPEDT* TAPEVL 

00900  COMMON/SVNM/  IPCF,IPRF 

01000  C  -  -  -  -  -  --------------------------- 

01100  C  WRITE  (2,2/  FLOW  ARRAYS  ON  DISK 
01200  IPNT=NX+(NZ-1 )*NXTL 

01300  CALL  ARRYMV ( 1,2,2) 

01400  CALL  STORIT (T APEPF , TPROF ,  I PRF ,  I PNT , 0 ) 

01500  INC=3-2*IFL0W 

01600  IGO~I CASE ( NZ ) 

01700  GO  TO  (100,200,300,300,400,500),  IGO 

01800  C  SET  POINTERS  FOR  NX=NXSTRT  EXCEPT  AT  END  OF  INITIAL  LINE  SWEEP 
01900  100  IPNTG=NX+(NZ+INC— 1 >*NXTL 

02000  GO  TO  650 

02100  C  SET  POINTERS  FOR  NX=NXSTRT  AND  END  OF  INITIAL  LINE  SWEEP 

02200  200  IF ( I FLOW. EQ. 2 )  GO  TO  210 

02300  I  PNT =NX+  (  NZ  1  —  1  )  -*NXTL 

02400  I PNT 1 1 =NX+NZ 1 *NXTL 

02500  IPNTG=NX+1+(NZ1-1 >*NXTL 

02600  GO  TO  220 

02700  210  IPNT=NX+ (NZ2-1 )*NXTL 

02800  I PNT 1 1=NX+(NZ2-2)*NXTL 

02900  i  riv7G=*NX-^i-*-iNZ2—  i  /*NXTL 

03000  220  CALL  STOR I T ( TAPEPF , TPROF , I PRF , I PNT , 1 ) 

03100  CALL  ARRYMV (2, 1 , 2) 

03200  CALL  ARRYMV (5, 1 , 2) 

03300  GO  TO  600 

03400  C  SET  POINTERS  FOR  GENERAL  CASE  AND  BEGINNING  OF  LINE  SWEEP 
03500  300  CALL  ARRYMV ( 6 , 1,2) 

03600  I PNTG=NX+ ( NZ+ I NC— 1 >*NXTL 

03700  NZNEXT=NZ+INC 

03800  IF( I  CASE ( NZNEXT ) . NE. 4 >  GO  TO  650 

03900  CALL  ARRYMV (6, 2, 1 ) 

04000  IPNT11*NX— 1+<NZ+2*INC— 1 >*NXTL 

04100  GO  TO  600 

04200  C  GENERAL  CASE... END  OF  LINE  SWEEP 

04300  400  IF(NX.EO.NXT)  RETURN 

04400  IF ( I FLOW. EQ. 2 )  GO  TO  410 

04500  I PNT =NX+ ( NZ 1  —  1  ) *NXTL 

04600  I PNT 1 1 =NX+NZ 1 *NX  TL 

04700  IPNTG=NX+1+(NZ 1— 1 )*NXTL 

04800  60  TO  420 

04900  410  IPNT=NX+(NZ2-1 >*NXTL 

05000  I PNT 1 1 =NX+ ( NZ2— 2 ) *NXTL 

05100  IPNTG=NX+1+ (NZ2-1  )*NXTl_ 

05200  420  CALL  STOR I T ( TAPEPF , TPROF , I PRF , I PNT , 1 ) 

05300  CALL  ARRYMV (2, 1,2) 

05400  CALL  ARRYMV (5, 1 , 2) 

05500  GO  TO  600 

05600  C  END  OF  ATTACHED  REGION 

05700  500  NZN=NZ 

05800  510  NZN=NZN+INC 

05900  NZNP=NZN+ INC 

06000  IF ( NSEP (NZNP).EO.l)  GO  TO  510 
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06100  IPNT=NX-1+NZN*NXTL 

06200  CALL  STOR I T  <  TAPEPF i TPROF i  I PRF i  I PNT » 1 ) 

06300  CALL  ARRYMV(2» li2) 

06400  CALL  ARRYMV (5i 1 i 2 ) 

06500  I  PNT  1 1  =NX—  l  +  (NZN-*-INC)  *NXTL 

06600  I PNTG=NX+NZN*NXTL 

06700  C  READ  IN  (1»1>  FLOW  ARRAYS  AND  <2. 2)  GEOMETRY  PARAMETERS 
06800  600  CALL  STOR IT  <  TAPEPF i  TPROF i I PRF i I PNT1 1 1 1 ) 

06900  CALL  ARRYMV ( 2» 1 » 1 ) 

07000  650  CALL  STOR I T ( TAPEGPi TPCFi I PCF i I PNTG i 1 ) 

07100  CALL  ARRYMV (3*2*2) 

07200  RETURN 

07300  END 


00 100  SUBROUTINE  SHIFT2 < NZ 1 , NZ2 > 

00200  INTEGER  TAPE  IN,  TAPEOT »  TAPEGP,  TAPEPF,  TAPEDT ,  TAPEVL 

00300  COMMON/BLC0/  NX iNZi NP»  NXT  » NZT ,  NTR,  NPT 1NXTL1 NZTL , NXSTRT i 

00400  *  NZSTRT,KC> IT, IFLOW, ICHORD, ISPAN, INTDIR 

00500  C0MM0N/BLC5/  F< 101 , 2, 2 ) , U < 1 01 , 2, 2 > , V < 101 > , G ( 101 , 2, 2 > , W< 101 , 2, 2 > , 

00600  *  T< 101 ) , TPROF ( 620 ) , TPCF  < 10 ) 

00700  COMMON/BLC9/  NSEP ( 4 1 ) , I  CASE ( 4 1 > 

00800  COMMON/TAPE/  TAPE  IN, TAPEOT, TAPEGP, TAPEPF, TAPEDT, TAPEVL 

00900  COMMON/SVNM/  IPCF,IPRF 

01000  C  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  — 

01100  IF ( I FLOW. EQ. 2 )  GO  TO  50 

01200  C  SEPARATION  ON  LOWER  BOUNDARY. .. IFLOW=l 

01300  IF ( NZ . GT . NZ 1 )  GO  TO  10 

01400  5  NZ 1 =NZ 1+1 

01500  NSEP ( NZ 1  )  =  1 

01600  IF ( NSEP (NZ1+1). EQ . 1 )  GO  TO  5 

01700  I PNT=NX— 1  +  < NZ 1  —  1 >*NXTL 

01800  CALL  STOR I T ( TAPEPF , TPROF , I PRF , I PNT , 1 ) 

01900  CALL  ARRYMV(2, 1,2) 

02000  I PNT 1 1 —NX— 1 +NZ 1 *NXTL 

02100  I PNTG=NX+ ( NZ 1  —  1 )*NXTL 

02200  I CASE ( NZ 1 ) =3 

02300  GO  TO  80 

02400  C  SEPARAR I ON  ON  UPPER  BOUNDARY. .. IFLOW=l 
02500  10  IF ( NZ . LT. NZ2 )  GO  TO  20 

02600  15  NZ2=NZ2- 1 

02700  NSEP  <  NZ2 ) - 1 

02600  IF ( NSEP ( NZ2— 1 ) . EQ. 1 )  GO  TO  15 

02900  Tcvwy. EO.NXT)  RETURN 

03000  I PNT =NX+ <  NZ 1 - 1 ) *NXTL 

03100  CALL  STOR I T  <  T APEPF , TPROF ,  I PRF,  I PNT , 1 ) 

03200  CALL  ARRYMV (2, 1,2) 

03300  IPNT1 l=NX+NZi*NXTL 

03400  IPNTG=NX+1+(NZ 1-1 )*NXTL 

03500  GO  TO  80 

03600  C  SEPARATION  AT  INTERIOR  POINT. .. IFLOW=l 
03700  20  NZN=NZ+ 1 

03800  NSEP  (  NZN-  1  )  ==  1 

03900  IF ( NSEP ( NZN+ 1 ) . EO. 1 )  GO  TO  20 

04000  IPNT=NX-1+<NZN-1 )*NXTL 

04100  CALL  STORIT  <  TAPEPF , TPROF , I PRF ,  I PNT , 1 ) 

04200  CALL  ARRYMV (2, 1,2) 

04300  IPNT1 1=NX— 1+NZN*NXTL 

04400  IPNTG-NX-MNZN-1  )*NXTL 

04500  I  CASE (NZN) =3 

04600  GO  TO  B0 

04700  C  SEPARAR I ON  ON  UPPER  BOUNDARY _ IFL0W=2 

04800  50  IF(NZ. LT. NZ2)  GO  TO  60 

04900  55  NZ2=NZ2— 1 

05000  NSEP ( NZ2 ) = 1 

05100  I F ( NSEP ( N Z 2— 1 ) . EO . 1 )  GO  TO  55 

05200  IPNT=NX  — 1-MNZ2-1  )*NXTL 

05300  CALL  STORIT( TAPEPF, TPROF, IPRF, IPNT, 1 ) 

05400  CALL  ARRYMV ( 2, 1 ,2) 

05500  I PNT 1 1 =NX— 1 + ( NZ2— 2 ) *NXTL 

05600  IPNTG=NX4-(NZ2-1  )*NXTL 

05700  I  CASE ( NZ2 ) — 3 

05800  GO  TO  80 

05900  C  SEPARATION  ON  LOWER  BOUNDARY. . . IFL0W=2 
06000  60  IF < NZ . GT . NZ 1 )  GO  TO  70 
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65  NZ1=NZ1+1 

NSEP ( NZ 1 >=1 

IF ( NSEP  <  NZ 1  +  1 ) *  EO. 1 )  GO  TO  65 
IF  <  NX • EQ • NXT )  RETURN 
IPNT=NX+<NZ2-1 >*NXTL 

CALL  STORI T ( TAPEPF i TPROF i IPRFi IPNT* 1 ) 
CALL  ARRYMV<2» 1 »2> 

IPNT1 1«NX+<NZ2-2)*NXTL 
IPNTG=NX+1+<NZ2-1 )«NXTL 
GO  TO  80 

C  SEPARATION  AT  INTERIOR  POINT. .  . IFL0W=2 
70  NZN=NZ-1 

NSEP ( NZN+ 1 )=1 

IF ( NSEP ( NZN— 1 ) . EQ. 1 )  GO  TO  70 
IPNT=NX-1+<NZN-1 >*NXTL 

CALL  STORITC TAPEPF i TPROF i IPRFi IPNT* 1 ) 
CALL  ARRYMV <2* 1*2) 

I PNT 1 1 -  NX— 1 +  <  NZN-2 ) *NXTL 
IPNTG=NX+(NZN-1 )*NXTL 
I  CASE ( NZN ) -3 
80  CALL  ARRYMV ( 5  *  1*2) 

CALL  STOR I T ( TAPEPF i TPROFi IPRF, IPNT11, 1 > 
CALL  ARRYMV ( 2i 1 i 1 ) 

CALL  STOR I T ( TAPEGP »  TPCF  i  IPCFi  IPNTGi  1 ) 
CALL  ARRYMV < 3f 2i 2) 

RETURN 

END 
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SUBROUTINE  SOLVEW 

INTEGER  TAPE  IN. TAPEOT  t  TAPEGP* TAPEPF,  TAPEDT » TAPEVL 
DIMENSION  ESAVE  < 101)i WSAVE < 101 ) 

COMMON /BLC0/  NXi NZi NP»NXT  .NZT.NTR*  NPT . NXTL  »NZTLi NXSTRT  i 

*  NZSTRT *  KC»  ITi  IFLOW,  I  CHORD.  ISPAN.  INTDIR 

COMMON /BLC3/  X (81 ) t  Z(41 >  « DETA ( 101  >  . ETA ( 101 ) » ETAE, VGP, CEL »  BEL  1  *  BEL 
COMMON/BLC5 /  F(101*2«2)*U(101*2*2>*V(101>*G(101*2*2>*W(101*2*2>* 

*  T  < 101 ) * TPROF ( 620 ) * TPCF ( 10) 

COMMON/ PRES/  PI (2i2)t P2 ( 2»  2 ) » P3 ( 2i 2  >  t P4 (2i2)i P5 (2i2) i 

*  P6(2f 2)»P7<2i  2)»P8(2»2)f  P9(2»2)f  P10(2»2) i Pll (2»2) t 

*  P12(2»2)iP13<2*2) i P14(2i2) 

COMMON/ PRSC/  PS 1 (2i2) iHl (2»2) » H2  <  2 »  2 ) »  UE ( 2 »  2 ) f 

*  WE (2*2) t CK1 (2 1 2 ) » CK2(2«2> «  CK1 2(2*2) i CK21 (2*2) * 

*  THETA (2*2) *  PFRS.UFRS. CNUFRS * UREF 1 . WNP 
COMMON/TURB/  E( 101 i 2i 2) i WT( 101 ) i WT2( 101 «  2*  2>  « ELT ( 101 ) 

COMMON/STOR/  CON* EPS* ICONVE* ITMAX »  WTEDG 

COMMON/ CONV/  ICOUNE*  ICOUNWt  INEGEi INEGW*  DELV* DELT 

COMMON/ CONS/  SALFA . SALFAS * SBET A >  SBETAS . SLAM DA . SS I GMA .  SSGMAS . USTOP 

*  Z I OT AE i  Z I OT AL  t  RW2 
COMMON /BLCE/  EDV<101) 

COMMON /BLC  7  /  PU  <  101  )  ,  PW  (  101  )  .  OU  <  101)*  O.U  (  101  ) 

COMMON /BLCY /  Y1 ( 101 ) 5  Y3 ( 101 ) i W1 i W2* W3* W4 
COMMON/ I PRT /  IPZi  IPX.  I  PRINT . EPSV. EPST 
COMMON /SHRT /  I SHORT 

COMMON/TAPE/  TAPE I N , TAPEOT »  TAPEGP i TAPEPF. TAPEDT . TAPEVL 
COMMON / REL X /  RFTRB. RFVEL 

DO  1  J =2 1 NP 

IF  v  U ( J i 2  »  2 ) • GE. . 999)  GO  TO  2 

1  CONTINUE 

2  NPL~ J 

NPM=NPL— 1 
U2=U ( NPL i 2 i 2 ) 

Ul»UCNPM«2*2) 

ET ADEL “ETA ( NPM ) +DET A  <  NPM ) *  < . 999-U1 >/<U2-Ul ) 

URUE=UREF 1 /UE ( 2 » 2 ) 

URUECT  =  2, *UREF 1 /UE ( 2* 2 ) *COS ( THETA ( 2* 2 )  ) 

REY=UE ( 2* 2 ) *PS 1 (2» 2) /CNUFRS 
UTAU4»V( 1 )**2 

IF  (  KC -  NE •  1  )  UTAU4=UTAU4+<URUE*T(  1  )  > **2+URUECT*V (  1 >*T( 1 > 

WTEDG-SORT ( REY*Z IOTAE/SALFAS ) / ( Z IOTAL*ETADEL ) 

UH AT  = ( UT AU4/ REY )  *■* . 25 
JSTOP— 1 

WCOEF-20 • /SBET A 
DO  25  J— 2 i NP 
UTOT»U( J*  2*  2) 

IF ( KC - NE .  1 )  UTOT —SORT ( UTOT **2+ ( URUE*W ( J*  2*  2 ) >**2 

*  +URUECT*UTOT*W( J*2* 2) ) 

UPLUS=UTOT /UHAT 

I F ( UPLUS • GT . USTOP )  GO  TO  26 
WT  <  J ) -WCOEF/ETA ( J )**2 
WT2  <  J  *  2*  2 ) =WT ( J ) **2 

25  JSTOP® JSTOP+1 

26  WT  < 1 ) ®WT  <  2 ) 

WT2  < 1 i 2* 2>»WT<2>**2 
E(lf2*2) ®0. 

I CONVE-0 
PU< 1 )=0. 

QU ( 1 ) =0 . 

PW< JSTOP )=WT2< JSTOP i 2.2) 
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QW< JSTOP>“0. 

E4M=1.+.5*SSGMAS*<EDV< 1 >+EDV<2> ) 

E5M* 1 . +• 5*SSIGMA*<EDV< 1 >+EDV<2> ) 

DO  150  J =2 1 NPM 
ESAVE ( J ) =E ( J i 2 »  2 ) 

WSAVE ( J ) =WT 2 <  J  >  2  *  2  ) 

E4P=1 .  + .  5*SSGMAS*  (  EDV  ( J  )  -+-EDV  <  J+ 1  )  ) 
E5P=1.+.5*SSIGMA*(EDV<J)+EDV( J+ 1 ) ) 

VB=V( J) 

TB=T ( J ) 

FB»F(J«2«2) 

UB=U< Jf 2f 2) 

6B- G ( J i 2»2) 

WB=W< J*2»2> 

EB=E< J, 2i 2) 

WTB=WT< J) 

E2-P1 <2, 2)*FB+P6<2*2)*GB 
DFB=CEL*(F< J»2i 2) -F(Ji  If  2) ) 

DGB=BEL 1 *  <  6  <  J i 2 1 2 ) -G < J i 2 1  1 )  ) 4 BEL2* < G < J i  1 1 1 )-6( Ji 1 1 2)  ) 
DVB=  <  E2+DFB+DGB ) / ( DET A(J)+DETA( J— 1 )  ) 
RETB=URUE**2*REY*EB/WTB 
GAMS= 1 . -COf\!*EXP(-RETB> 

SHEAR2«VB*VB 

IF  <  KC. NE.  1 )  SHEAR2-SHEAR2+ ( URUE*TB ) *-*2+URUECT*TB*VB 

PR0DE=GAMS*REY*SHEAR2/WTB 

DISSE=SBETAS*WTB 

ABE*.  3 

PS I F I X=PRODE/D I SSE— .  7 
IMKblMX.bl  .  ABt;  Abt  =  Pbih  i  a 
X3=PR0DE—  (  1 . 4-ABE  )  *D  I  SSE 

X7=-(CEL*UB+BEL2*WB)*E< J» 1,2) —BEL 1 *WB*E (  J »  2  ? 1 > 

*  +BEL2*WB*E< J, 1,1) 

A 1 =  E4M*Y3 ( J ) —DVB 

B 1 =—E4P*Y 1 ( J ) -E4M  w  <  J ) +X3- ( CEL*UB+BEL 1 *WB ) 

Cl=  E4P*Y1  < J)+DVB 
D1 *X7- ABE*D I SSE*EB 
B 1 S=B 1 + A  1 *QU ( J— 1 ) 

D1S=D1--A1*PU< J-l ) 

PU( J)*DlS/eiS 

QU ( J ) =—C 1 /BIS 

IF< J.LE. JSTOP)  GO  TO  145 

WT2B*WT2< Ji 2 *  2  > 

GAM=SALFA*< 1 . -CON*EXP ( -RETB/RW2 > ) /GAMS 
DLDY==(ELT(  J+l  )-ELT(  J-l  )  )  /  (  DET  A  (  J  )  +DET  A  (  J-  1  )  > 
DLDY2=DLDY*  DLDY 
PR0DW=GAM*PR0DE+X9 

SBET AT “SBET A+2  •  *SS  I  GMA*REY*DLDY2-*URUE*-*2 

DISSW=SBETAT*WTB 

ABW*.  3 

PS IF I X*PRODW/D I SSW~ .  7 
IF( PS  IF  I  X. GT. ABW)  ABW*PSIFIX 
PSISTP*. 5*PR0DW/D I SSW— . 25 
IF <  PS I ST P. GT. ABW )  ABW^PSISTP 
X6*PR0DW- ( 1 . +ABW  >*DISSW 

X8=-  <  CEL*UB+BEL2-*WB  )  *  WT2  (  J  i  1 1  2  )  -BEL  1  *WB*WT2  (  J*  2,  1  ) 

*  +BEL2*WB*WT2 (  J ,  1 1  1  > 
X9=2.  *<P13<2,2>  *UB+P14  <  2,2)  *WB  ) 

A2»  E5M*Y3( J)-DVB 

B2=-E5P*Y 1 ( J ) — E5M* Y3 <  J  >  +X6- <  CEL*UB+BEL 1 * WB ) 

02*  E5P*Y1  <  J)-*-DVB 
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12100 
12200 
12300 
12400 
12500 
12600 
12700 
1 2800 
12900 
13000 
13100 
13200 
13300 
13400 
13500 
13600 
13700 
13800 
13900 
14000 
14100 
14200 
14300 
14400 
14500 
14600 
14700 
14800 

15000 

15100 

15200 

15300 

15400 

15500 

15600 

15700 

15800 

15900 

16000 

16100 

16200 

16300 

16400 

16500 

16600 

16700 

16800 

16900 

17000 

17100 

17200 

17300 

17400 


D2=X8-ABU*DISSU*WT2B 
B2S=B2+A2*QW< J-l ) 

D2S=D2-A2*PW<  J-l ) 

PW< J)=D2S/B2S 
QW ( J ) =— C2/B2S 
E4M=E4P 
E5M=E5P 

WT2EDG=WTEDG*WTEDG 
E  <  NPL i2i2)“ZI OT AE 
WT2(NPL, 2i 2>=WT2EDG 
KON— NPM 

DO  170  J NPM 

E  <  KON i2i2)=PU( KON ) +QU ( KON ) *E ( KON+ 1,2.2) 

IF (KON. LE. JSTOP)  GO  TO  170 

WT2 ( KON , 2 *  2 ) =  PW ( KON ) +GW ( KON ) *WT2 ( KON+ 1,2,2) 

KON-KON-1 
I COUNE— 0 
INEGE=0 

DO  175  J=2 , NPM 

IF(E( J, 2, 2) . GE. 0. )  GO  TO  176 

I NEGE= I NE6E+ 1 

E( J, 2, 2)=. 5*ESAVE< J) 

I F ( J . GE •  ( NPL— 3 )  )  E  (  J ,  2 , 2  >  =  Z  I OT AE 
GO  TO  175 

I F ( ABS ( E  <  J , 2 , 2 ) -ESAVE ( J  >  )  .  GT. EPS*ESAVE ( J )  )  I COUNE** I COUNE+ 1 
E( J, 2, 2 ) =ESAVE ( J ) +RFTRB* ( E ( J, 2, 2 ) —ESAVE ( J ) ) 

I COUNW=0 
INEGW- 0 

DU  1V!D  <J— «JS*I  OP,  NPM 

IF(WT2( J, 2, 2) . GT. 0. )  GO  TO  196 

WT2( J, 2, 2)-. 5*WSAVE< J) 

I NEGW= I NEGW+ 1 
GO  TO  197 

I F ( ABS ( WT2 ( J , 2 , 2 ) - W5AVE  <  J )  ) . GT . EPS*WSAVE ( J )  )  I COUNW= I COUNW+ 1 
WT2 ( J, 2, 2) =WSA VE ( J ) +RFTRB*  < WT2 <  J »  2 ,  2  ) -WSAVE ( J )  ) 

WT( J>=SQRT(WT2< J, 2,2) ) 

CONTINUE 

I  COUNT  =  I  COUNU+  I  COUNE  +  I  ,.EGW+  I NE GE 
I F ( I COUNT. LE. 2 )  ICONVE=l 
NPP=NPL+ 1 
DO  200  J=NPP , MPT 
E( J, 2, 2 ) =E ( NPL , 2,2) 

WT ( J ) =WTEDG 

WT2 ( J  f  2 1  2 ) =WT2EDG 

IF  < I SHORT • EO.  1  )  RETURN 

I F ( I  PR  I NT • EO . 0 )  RETURN 

IF ( IT. EO . 1 )  WRITE (TAPEOT, 40) 

WRITE (TAPEOT, 50)  ITiV(l), DELV, T<  1  ) , DELT,  I NEGE  * 

*■  INEGW, I COUNE, ICOUNW 

RETURN 

FORMAT ( 1H0, 1X2HIT, 5X7HVCWALL ) , 7X7HDELV< 1  * . 7*  ^  '  , 

*■  7X5H  I  NEGE ,  3X5HINEGW*  3*6HI  20UNE  ,  .  < I  ' 

FORMAT (1H  , I 3 , 1 P4E 1 4 . 5 , 4 I 8 ) 

END 
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MICROCOPY  RESOLUTION  TEST  CHART 
NATIONAL  BUREAU  OF  STANDARDS-  1963-A 


00100  SUBROUTINE  S0LV6 

00200  DIMENSION  USAVE ( 101 ) . WSAVE ( 101 ) 

00300  COMMON/BLC0/  NX » NZ .  NP. NXT »NZT» NTR. NPT i NXTL»  NZTL i NXSTRT . 

00400  *  NZSTRT.KC. IT. IFLOW. I  CHORD. ISPAN. INTDIR 

00300  C0MM0N/BLC3/  X(81)»Z<41)»  DETA( 101  > . ETA < 101 > » ETAE. VGP. CEL » BEL 1  *  BEL 

00600  C0MM0N/BLC5/  F( 101 .2*2) » U< 101,2.2) .V( 101 )i6( 101.2,2),W< 101.2.2) , 

00700  *  T( 101 ) i TPROF  <  620 ) »  TPCF ( 10 ) 

00800  C0MM0N/BLC7/  PU< 101 ) , PW< 101 ) , OU< 101 ) »QW< 101 ) 

00900  COMMON/BLCY/  Y1 ( 101 > . Y3< 101 ) . W1 . W2. W3. W4 

01000  COMMON/ I PRT /  IPZ. IPX. I PR I NT . EPSV. EPST 

01100  COMMON/ CONV /  I COUNE . I COUNW . I NEGE . I NEGW i DEL V . DELT 

01200  COMMON/RELX/  RFTRB. RFVEL 

01300  C  —  — - - - - - - - - - 

01400  VWSAVE=V ( 1 ) 

01500  TWSAVE=T(1> 

01600  NPM=NP— 1 

01700  KON*NPM 

01800  DO  10  J-2.NPM 

01900  USAVE ( KON ) =U ( KON .2.2) 

02000  WSAVE  <  KON ) =W  <  KON .2.2) 

02100  U ( KON , 2 , 2  >  =PU (  KON )  +QU  ( KON >  *U  (  KON+ 1.2.2) 

02200  U  (  KON .  2  .  2 )  =PW  ( KON )  +«3W  ( KON )  »W  ( KON+ 1.2.2) 

02300  10  KON=*KON— 1 

02400  U(1.2.2) =0 . 

02500  W<1,2,2>=0. 

02600  DO  20  J=2»NPM 

02700  U ( J . 2 . 2 ) =US A VE ( J ) +RFVEL* ( U ( J  »  2 . 2 ) -USAVE  <  J ) ) 

02800  20  W<J.2.2) =WSAVE  <  J )+ RFVEL*  <W(J»2»2) —WSAVE  <  J ) ) 

02V00  F< 1.2.2 )“»- 

03000  G ( 1 . 2. 2 ) =0. 

03100  CF1=0. 

03200  CG 1 =0 . 

03300  DO  30  J=2,NP 

03400  CF=U(J,2,2) 

03500  CG=W ( J. 2. 2 ) 

03600  F(J»2,2)=F(J-1,2,2>+.5*(CF+CF1 )*DETA(J-1 ) 

03700  G( J. 2. 2>=G( J-l .2.2)+. 5*( CG+CG1 )*DETA< J— 1 ) 

03800  CF1=CF 

03900  30  CG1=CG 

04000  DO  35  J=2»  NPM 

04100  V<  J)s=(U(  J+l  .2.  2)— U(  J— 1 .2.2)  ) / ( DETA ( J ) +DETA (  J-l  )  ) 

04200  35  T<J)=(W(J+1,2,2>-W( J-1,2,2) )/( DETA ( J > +DETA ( J-l ) ) 

04300  V(NP ) ( 1 .  — U ( NPM.  2.2)  ) /DETA  (NPM) 

04400  T ( NP )=(W(NP,2,2) — W ( NPM. 2.2) ) /DETA (NPM) 

04500  V( 1 )«— W1*U( 1.2.2 )+W2*U( 2.2.2 >-W3*U( 3.2.2 )+W4*U(4, 2. 2) 

04600  T ( 1 >=-Wl*W( 1,2,2)+W2*W(2,2,2)-W3*W(3,2,2)+W4*W(4,2,2) 

04700  DELV*V ( 1 ) — VWSAVE 

04800  DELT=»T  ( 1 )  -TWSAVE 

04900  RETURN 

05000  END 


00100  subroutine  start 

00200  INTEGER  TAPE IN*  TAPEOT*  TAPEGP*  TAPEPF t  TAPEDT »  TAPEVL 

00300  DIMENSION  VPLUS(101) 

00400  COMMON /BLC0/  NX* NZ . NP. NXT. NZT. NTR* NPT. NXTL. NZTL. NXSTRT. 

00900  •  NZSTRT.KC. IT. IPLOW. ICHORD* ISPAN* INTDIR 

00600  C0MM0N/8LC3/  X (01 ) * Z <41 ) . DETAC 101 ) » ETA( 101 > . ETAE. VGP* CEL. BEL1 » BEL 

00700  C0MH0N/BLC3/  F< 101 . 2. 2> . U< 101 » 2* 2 ) . V< 101 ) , G< 101 . 2. 2) . W< 101 » 2* 2) * 

00800  *  T  < 101 ) »  TPROF ( 620 ) »  TPCF ( 10 ) 

00900  COMMON/ PRSC/  PS1 < 2* 2) » HI < 2* 2) * H2( 2. 2) . UE( 2* 2> . 

01000  •  UE(2*  2) «  CK1 (2*  2) *  CK2(2*  2) *  CK12<2*  2 >  *  CK21 (2*2)* 

01100  *  THETA <2*  2) . PPRS*  UFRS*  CNUFRS. UREF1 , WNP 

01200  COMMON/BLCE/  EDV(101> 

01300  COMMON/TURB/  E< 101  *  2* 2 ) * WT < 101 ) » WT2< 101 . 2. 2 ) . ELT( 101 ) 

01400  COMMON/STOR/  CON. EPS* I CONVE. I TMAX* WTEDG 

01900  COMMON/TAPE/  TAPEIN* TAPEOT* TAPEGP. TAPEPF* TAPEDT. TAPEVL 

0 1 600  COMMON/ CONS/  SALF A *  SALFAS  *  SBETA • SBETAS . SLAMDA  *  SS I GM A «  SSGMAS . USTOP 

01700  *  ZIOTAE* Z TOTAL. RW2 

01800  COMMON/ I NTG/  ALAMI <41 ) . CFXI (41 ) . CFZI <41 ) , DELTAI (41 ) . 

01900  *  THETXI (41 ) .THETZI <41 > 

02000  DATA  PI/3. 14199263/* AKAPPA/. 41/ 

02100  C - - - - - 

02200  PS1 ( 1.21-PS1 <2. 2) -HI (2»2)*<X(NXSTRT+1 >-X< NXSTRT) ) 

02300  ALAMM-ALAM I ( NZ ) 

02400  CFX-CFXI (NZ) 

02300  CFZ-CFZI(NZ) 

02600  DELTA»DELT  A I ( NZ ) 

02700  THETAX«THETX I ( NZ ) 

02800  THETAZ-THETZ I ( NZ  > 

S29SS  RES  »OEi  2*  2 1  l  a*/ J  /CNUt-KS 

03000  REDELT-UE (2*2) * DELTA/ CNUFRS 

03100  ROOTRS-SORT(RES) 

03200  URUE“UREF 1 /UE (2*2) 

03300  URUE2»URUE*«2 

03400  UTAU2».9*(CFX+ABS(CFZ) ) 

03500  UTAU-SORT ( UTAU2 > 

03600  DELTAP-UREF 1 *UTAU*DELTA/ CNUFRS 

03700  YSCALP«UTAU*URUE*ROOTRS 

03800  C0STH»C0S(THETA<2.2) ) 

03900  SINTH-SIN(THETA<2.2> > 

04000  COTTH-COSTH/SINTH 

04100  DO  10  J»1.NP 

04200  10  YPLUS< J)-YSCALP*ETA(J> 

04300  C  VELOCITY  PROFILES 
04400  D0TXM3=DELTA/THETAX-3. 

04500  D0TZM3=3. 

04600  IF( ABS ( THETAZ ) . GT. 1 . E— 8 )  D0TZM3=DELTA/THETAZ-3. 

04700  ENX».5*(D0TXM3+S0RT(D0TXM3*D0TXM3-8. )) 

04800  ENZ*.  5*  (  D0TZM3-+ SORT  ( D0TZM3*D0TZM3— 8.  )  ) 

04900  SUBSCL».5*R00TRS*URUE2 

05000  OTSCLX* ( 1 . +WE (2.2) *COSTH/UE ( 2. 2 ) ) * ( ROOTRS/REDELT ) **  < 1 . /ENX > 

05100  OTSCLZ*  < WE  <  2 . 2 ) *SINTH/UE ( 2. 2 ) ) * ( ROOTRS/REDELT ) **  < 1 .  /ENZ ) 

05200  V< 1 >»SUBSCL^CFX 

05300  T< 1 >«SUBSCL*CFZ 

05400  U( 1 1 2*  2 ) =0. 

05500  W( 1  *  2*  2  >  *0. 

05600  DO  20  J»2*NP 

05700  I F ( YPLUS ( J ) . GE . DELT AP )  GO  TO  30 

03800  U( J«  2*  2 ) -CFX*SUBSCL*ETA  <  J ) 

03900  W( J*2«2)»CFZ*SUBSCL*ETA( J) 

06000  UOUT-OTSCLX*ETA < J)##< 1 . /ENX  > 


06100 
06200 
06300  20 

06400  30 

06500 
06600 
06700 
06800 
06900 
07000  35 

07100 
07200  40 

07300 
07400 
07500  C 
07600 
07700 
07800 
07900 
08000 
08100 
08200 
08300 
08400 
08500 
08600  50 

08700  C 
08800 
08VW0 
09000 
09100 
09200 
09300 
09400  65 

09500 
09600 
09700  68 

09800  70 

09900  C 
10000 
10100 
10200 
10300 
10400 
10500 
10600 
10700 
10800  75 

10900 
11000 
11100 
11200 

11300  78 

11400  80 

11500 
11600 
11700 
11800 
11900 


WOUT«OTSCLZ*ETA  <  J ) **  < 1 . /ENZ ) 

IF<U<Ji2f2>.GT. UOUT )  U  <  J  t  2 1  2  >  «*UOUT 
IF<ABS<W<J1 2*2) ) . GT • ABS  <  WOUT ) )  W< J* 2t 2 ) -WOUT 
WSCALE* 1 . / <  URUE*S I NTH ) 

DO  40  J»2iNP 

IF <  YPLUS  <  J ) • GE.  DELTAP )  GO  TO  35 
U< Jf 2*2>»U< Jf 2f 2>-W< Jf 2t2>*C0TTH 
W< J»2i 2>»WSCALE*W< J* 2t 2) 

GO  TO  40 
U<J*2f2)«l. 

W < J  i  2 1 2 ) ®WE < 2 1 2 ) /U REF 1 
CONTINUE 

V<1 >«V< 1 >— T < 1 )*COTTH 
T  < 1 ) ~WSCALE*T  < 1 ) 

STREAMFUNCTION  COMPONENTS 
F  < 1 »  2t  2 ) —0. 

G  <  1  ?  2*  2 ) =0. 

CF1«0. 

CGI*5©. 

DO  50  J=2»  NP 
CF«U< Ji2i2) 

CG—W<  J i  2j  2 ) 

F< Jf  2i  2 ) =F  < J— 1 i2i2)  +  .5#( CF+CF1 )*DETA< J— 1 ) 
G<J*2i2)»G< J-lf 2i 2 >+.5* <CG+CG1 >*DETA< J-l > 

CF1=CF 

CG1«CG 

TURBULENT  MIXING  ENERGY 

ESCALE=ALAMM*UTAU2/SALFAS 
fcWAICH«ESCALE*<COS<5.*PI/DPI.  TAP)  ) **2 
DO  70  J»1,NP 

IF<YPLUS< J) .GT. 10. >  GO  TO  65 
E<Ji2f2) =EMAT  CH*  < •  1*YPLUS  <  J  ) >**4 
GO  TO  70 

IF < YPLUS <J).GE. DELTAP)  GO  TO  68 

E <  Jf  2i 2 >  =ESCALE* <  COS < . 5*PI *YPLUS <  J ) /DELTAP )  >  **2 
IF<E< Ji2f 2>.GT. ZIOTAE)  GO  TO  70 
E<J*2*2) =  Z IOTAE 
CONTINUE 

TURBULENT  DISSIPATION  RATE 

YPLSUB=20. *SALFA5*AKAPPA/SBETA 
WSUB*=<20.  /SBETA )  **2 

WWAL=  <  ALAMM*YSCALP/ <  SALFAS*AKAPPA ) >  **2 
WWAK=  <  PS1 < 1 ,2>*URUE/<ZI0TAL*DELTA) >**2/SALFAS 
DO  80  J»2f NP 

IF < YPLUS <J).GT. YPLSUB)  GO  TO  75 
WT2 ( Jf  2i 2) ~WSUB/ETA  <  J ) **4 
GO  TO  80 

WTEST 1 *WWAL/ETA  <  J ) **2 
WTEST25*  WWAK*E  <  Jt  2*  2) 

IF < WTEST 1 . LT. WTEST2. OR. YPLUS < J ) • GE. DELTAP )  GO  TO  78 
WT2 <  Jv  2i  2 ) **WTEST  1 
GO  TO  80 

WT2< J,2f 2)«WTEST2 
WT  <  J ) *SQRT  <  WT2  <  J 1 2 1  2 ) ) 

WTEDG»WT<NP) 

WT  < 1 ) *WT  <  2 ) 

WT2< If 2»2)-WT<2)**2 

RETURN 

END 
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00100  SUBROUTINE  STORIT < I UN IT *  A*  NDIM* I PNT » N) 

00200  C  it*********************-************************************** 
00300  C  THIS  SUBROUTINE  COORDINATES  TEMPORARY  STORAGE  AND  RETREIVAL 
00400  C  OF  DATA  FROM  DISK.  DATA  ARE  STORED  AS  A  SERIES  OF  BLOCK 
00S00  C  DATA  IN  THE  A  ARRAY  WHOSE  DIMENSION  IS  NDIM.  DURING  A  READ 
00600  C  FROM  THE  SUBROUTINE  (N«l>»  THE  FIRST  (IPNT-1)  BLOCKS  ARE 
00700  C  SKIPPED  AND  THE  A  ARRAY  CORRESPONDING  TO  I PNT  IS  READ  FROM 
00800  C  DISK.  DATA  ARE  WRITTEN  ON  DISK  WHEN  N»0. 

00900  C  ************************************************************ 
01000  INTEGER  TAPEIN*  TAPEOT *  TAPEGP*  TAPEPF *  TAPEDT* TAPEVL 

01100  DIMENSION  A ( 620 ) 

01200  COMMON/TAPE/  TAPEIN* TAPEOT. TAPEGP. TAPEPF* TAPEDT. TAPEVL 

01300  C  ****************** 

01400  C  WRITE  DATA  ON  DISK 
01500  C  ****************** 

01600  C 

01700  IF (N. NE. 0 )  GO  TO  10 

01800  C 

01900  C*****IBM  OPTION***** 

02000  C  WR I TE  < I UN I T  *  I PNT )  A 

02100  C*****CDC  7600  OPTION***** 

02200  C  CALL  WR I TMS ( I UN I T »  A  »  ND I M i I PNT ) 

02300  C*****UNI VAC  1108  OPTION***** 

02400  INDEX=NDIM*( IPNT-1 ) 

02500  CALL  NTRAN(IUNIT. 10) 

02600  CALL  NTRANC I UNIT* 6. INDEX) 

02700  CALL  NTRAN ( I UN I T  *  1  *  ND I M  *  A  *  LFLAG ) 

02800  1  I F ( LFLAG . EQ . — 1 )  GO  TO  1 

(929(0 Id  IF  (LFLAG.  LT.  — 1  )  GO  TO  20 

03000  RETURN 

03100  C 

03200  C  ******************* 

03300  C  READ  DATA  FROM  DISK 
03400  C  ******************* 

03500  C 

03600  10  CONTINUE 

03700  C 

03800  C*****IBM  OPTION***** 

03900  C  READ  < I UN IT * I PNT)  A 

04000  C*****CDC  7600  OPTION***** 

04100  C  CALL  READMS ( I UN IT *  A*  NDIM* I PNT) 

04200  C*****UN I VAC  1108  OPTION***** 

04300  INDEX=NDIM*( IPNT-1 ) 

04400  CALL  NTRAN  < I UN I T * 1 0  > 

04500  CALL  NTRAN ( IUN IT. 6. INDEX) 

04600  CALL  NTRAN ( I UN IT *  2*  NDIM* A*  LFLAG) 

04700  11  I F  <  LFLAG .  E»i .  —  1  )  GO  TO  1 1 

04800  IF (LFLAG. LT. — 1 )  GO  TO  20 

04900  RETURN 

05000  C 

05100  20  IF (N. EG. 0 )  WRITE (TAPEOT *  5 ) 

05200  IF(N.EG.l)  WRITE (TAPEOT* 6) 

05300  IF < LFLAG. EG. -2)  WRITE (TAPEOT*  2)  LFLAG 

05400  IF(LFLAG. EG. -3)  WRITE (TAPEOT. 3)  LFLAG 

05500  IF (LFLAG. EG. —4)  WRITE (TAPEOT* 4)  LFl  AG 

05600  STOP 

05700  C  ********************************************************** 
05800  2  FORMAT ( 1H  »  1 X7HLFLAG  =•*  13. 5X23HEND  OF  FILE  ENCOUNTERED) 

05900  3  FORMAT ( 1H  » 1 X7HLFLAG  *  *  13*5X1 2HDE V I CE  ERROR) 

06000  4  FORMAT (1H  « 1X7HLFLAG  ■* 13* 5X20HTRANSMISSION  ABORTED) 
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REM  #########*##*##■»###*#######•»#*##################### 
REM  *******  INTERPOLATION  ROUTINE  TO  ACCOUNT  ******* 
REM  *******  FOR  EFFECTS  OF  ******* 

REM  *******  DISPLACEMENT  THICKNESS  ******* 

REM  *************************************************** 
DEFUSR 1 “&H7D00 
CLEAR  1000 

DIM  X(41)tZ(21)t DISP ( 41 v  21 ) 

DIM  ZB  <  41 1 21 ) 

DIM  XRC41 ) 9 YRC41 vll)fZR(ll>« YN<41 til) 

OPEN  "1*9  If "DISPL" 

INPUT  "NXTL  =  “;NX 
INPUT  "NZTL  «MNZ 
FOR  1*1  TO  NX 
FOR  J=1  TO  NZ 
INPUT#lt  DISP ( I i J ) 

NEXT  J 
NEXT  I 
CLOSE 

INPUT  "FILE  NAME  FOR  NONORTHOGONAL  MESH" ?F$ 

OPEN  "IM»F$ 

FOR  J=1  TO  NZ 
INPUT#lt  H$»Z(J) 

FOR  1=1  TO  NX 

INPUTttlt  D»  X< I) 1D1 ZB(If J)»DtDf Df Di DiDf D»D 
NEXT  I 
NEXT  J 
CLOSE 

OPEN  “  i  ::  t  i  t  c  HESHPRU 
INPUT  "IMAX  =  "?IM 
INPUT  "JMAX  =  "?JM 
FOR  1=1  TO  IM 
INPUT#li  XR( I ) 

NEXT  I 

FOR  J=1  TO  JM 
INPUT#lt  ZR ( J ) 

NEXT  J 

FOR  1=1  TO  IM 
FOR  J»1  TO  JM 
INPUT#1 i  YR(ItJ) 

NEXT  J 
NEXT  I 
CLOSE 

NAME  " MESHPR"  AS  "MESHPR/OLD" 

YN»USR1 (Xi ZtDISPt ZBt  XRt  YRt ZRi NX t NZ t IMt  JM ) 

OPEN  "O" i 1» "MESHPR" 

FOR  1=1  TO  IM 

PRINT#1  USING  "##.###### tttt-iXRCI) 

NEXT  I 

FOR  J=1  TO  JM 

PR I NT#1  USING  ■##. ######* ttt ■ , ZR( J) 

NEXT  J 

FOR  1=1  TO  IM 
FOR  J«1  TO  JM 

PRINT#1  USING  YN< It J) 

NEXT  J 
NEXT  I 
CLOSE 
END 
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